Importer les données

library(datasets)
data(iris)


Exploration de données

Avant de débuter des analyses, il est important de se familisariser avec son jeu de données afin d’avoir une idée de sa structure. Cette étape permet d’identifier des motifs, des tendances ou des relations préalablement aux tests statistiques, puis de vérifier la qualité des données. Elle permet aussi d’évaluer si les données ont besoin d’être transformées avant de procéder aux analyses.

Quelques fonctions de base

Pour vérifier le nombre de lignes et de colonnes de votre dataset:

dim(iris)
## [1] 150   5

Pour voir les 6 premières lignes du jeu de données:

head(iris)
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa

Pour connaître le nom des colonnes:

names(iris)
## [1] "Sepal.Length" "Sepal.Width"  "Petal.Length" "Petal.Width"  "Species"

Pour vérifier la classe de chaque variable dans la table de données:

str(iris)
## 'data.frame':    150 obs. of  5 variables:
##  $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
##  $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
##  $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
##  $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
##  $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...

Pour obtenir certaines statistiques descriptives de base, telles que le minimum, le maximum puis la moyenne (variables continues), ainsi que le nombre d’observations (variables catégoriques):

summary(iris)
##   Sepal.Length    Sepal.Width     Petal.Length    Petal.Width   
##  Min.   :4.300   Min.   :2.000   Min.   :1.000   Min.   :0.100  
##  1st Qu.:5.100   1st Qu.:2.800   1st Qu.:1.600   1st Qu.:0.300  
##  Median :5.800   Median :3.000   Median :4.350   Median :1.300  
##  Mean   :5.843   Mean   :3.057   Mean   :3.758   Mean   :1.199  
##  3rd Qu.:6.400   3rd Qu.:3.300   3rd Qu.:5.100   3rd Qu.:1.800  
##  Max.   :7.900   Max.   :4.400   Max.   :6.900   Max.   :2.500  
##        Species  
##  setosa    :50  
##  versicolor:50  
##  virginica :50  
##                 
##                 
## 

Pour connaître les niveaux d’une variable catégorique:

levels(iris$Species)
## [1] "setosa"     "versicolor" "virginica"


Visualisation d’un jeu de données

La fonction plot() est une fonction versatile qui permet de créer une grande variété de figures avec les données brutes. Dépendemment du type de donnée et des arguments fournis à la fonction, plot() peut générer plusieurs types de figures comme:

Une seule matrice de scatterplots figure peut être créée avec toutes les variables.

plot(iris)

Pour créer un scatterplot de la relation entre une variable spécifique et une autre, il suffit d’entrer celles-ci dans la fonction, séparées par le signe ~. Par exemple, si on veut visualiser la relation entre la longueur des sépales et celle des pétales dans le jeu de données iris, on peut l’écrire ainsi:

plot(iris$Sepal.Length~iris$Petal.Length)

À l’aide de la fonction boxplot(), on peut visualiser la dispersion d’une variable selon les groupes d’une variable catégorique. Ce type de figure permet aussi d’identifier rapidement des valeurs aberrantes ou des anomalies.

boxplot(iris$Petal.Length~iris$Species)

La fonction hist() génère un histogramme de la variable qui y est précisée.

hist(iris$Sepal.Width)


Conseils pour des figures et tableaux réussis

Les figures sont des représentations visuelles des résultats. Elles rendent la lectures des résultats principaux plus facile et permettent de mettre en évidence des tendances ou motifs intéressants. On veut pouvoir la comprendre sans avoir à osciller entre la figure et le texte. Chaque figure devrait correspondre à un message général de votre étude. Chacune devrait répondre en totalité ou en partie à un objectif.

Voici quelques aspects à considérer pour présenter vos figures:

  • Un titre descriptif: doit contenir les variables mesurées, les unités de mesure, le nom commun et en latin du taxon (si applicable). Le titre doit fournir assez d’information pour qu’on comprenne le contexte de la figure sans devoir se référer au texte.

  • Titre des axes: doivent comprendre les variables et leur unité de mesure.

  • Les données: les données brutes (non transformées!) doivent être présentées.

  • Barres d’erreurs, intervalles de confiance et/ou de prédiction: pour un graphe à barres, ajouter les barres d’erreur afin de représenter la variabilité associée aux données. Pour les figures de régression, il est important d’ajouter l’intervalle de confiance ou de prédiction (dans le cas des résultats d’un modèle linéaire généralisé binomial).

  • Légende dans la figure: une petite légende peut être nécessaire pour distinguer les traitements (couleur, type de ligne, etc.)

  • Indice de significativité: on doit retrouver la valeur p d’une relation, des astérisques au-dessus de graphes à barres pour indiquer si la relation est significative, le R2 d’un modèle.

Les tableaux sont un bon choix pour présenter de l’information numérique détaillée. Ils présentent habituellement des résultats plus complexes qui seraient trop encombrants à inclure dans une figure ou dans le texte. Généralement, si les données ne peuvent être présentées en une ou deux phrases, un tableau est nécessaire. Les lignes et les colonnes doivent contenir le nom de la variable ainsi que l’unité de mesure. Un tableau résumé des statistiques peut inclure, par exemple, la moyenne, l’écart-type, les intervalles de confiances, les degrés de liberté, la valeur p et autres statistiques (comme la F value).

Légende d’une figure ou d’un tableau

Les légendes servent à compléter l’information qui est présentée dans le tableau ou la figure. On y retrouve entre autre les tailles d’échantillons (n), valeur p, des descriptions des abréviations, la méthode de collection, le nombre de réplicats, etc.

Annexe ou pas?

Les données ou les figures qui ne contribuent pas directement à l’histoire principale de votre rapport peuvent être rassemblées en annexes. Dans cette section, on peut retrouver, par exemple, des cartes du site d’échantillonnage, du matériel utilisé pour récolter les données, des tableaux avec davantage d’informations sur les modèles (pas de capture d’écran du summary!), etc.


Coder des graphiques avec ggplot

Avant de commencer, identifiez vos variables dépendantes et indépendantes. Puis, déterminez quels sont les types de ces variables avec lesquelles vous allez travailler.

Sont-elles catégoriques (et nominales ou ordinales), ou bien numériques (et continues ou discrètes)?

C’est ce qui déterminera le choix de type de figure pour représenter vos données - un histogramme, un scatterplot, un boxplot, un barplot ? Autre?

Consultez ce lien pour mieux comprendre:

https://statisticsbyjim.com/basics/data-types/


ANOVA

La fonction de base plot() permet de visualiser rapidement un jeu de données, par exemple avec un histogramme. Par contre, vous verrez que son utilisation peut devenir limitée lorsqu’il s’agit de réaliser des figures plus complexes ou simplement modifier certains paramètres graphiques.

C’est pour cela que nous suggérons d’utiliser la fonction ggplot(), du package ggplot2. Cette fonction permet de réaliser des graphiques de façon plus intuitive et permet de les mettre en page plus facilement qu’avec la fonction plot(). Elle est aussi mieux documentée, il est donc plus facile de comprendre et utiliser les différents aspects de la fonction. Google est d’ailleurs votre meilleur allié pour la réalisation de vos figures!

Contrairement à la fonction plot(), ggplot() fonctionne par couches. Une figure ggplot() commence avec la fonction ggplot(). Elle sert à “préparer” la figure: on spécifie le jeu de données à utiliser, puis on choisit les variables qui formeront nos axes.

La fonction ggplot() nécessite deux arguments: le dataset (jeu de données), puis l’argument aes(). Ce dernier nous permet d’assigner des variables du dataset aux composantes du graphique (par exemple, les axes x et y). Voici un exemple, toujours avec le jeu de données Iris. Nous allons tester si la largeur des sépales diffère entre les espèces.

(Revoir l’atelier 3 pour l’explication de l’ANOVA, incluant les postulats et critères d’utilisation, que je passe ici)

library("ggplot2")

head(iris)
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa
boxplot(Sepal.Width ~ Species,
  data = iris)

On produit d’abord le modèle linéaire, et le test d’ANOVA, avec Sepal.Width comme variable dépendante et Species comme variable indépendante.

modele.SW<-lm(Sepal.Width ~ Species, data = iris)

#Test anova:
anova(modele.SW)
## Analysis of Variance Table
## 
## Response: Sepal.Width
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## Species     2 11.345  5.6725   49.16 < 2.2e-16 ***
## Residuals 147 16.962  0.1154                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(modele.SW)
## 
## Call:
## lm(formula = Sepal.Width ~ Species, data = iris)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.128 -0.228  0.026  0.226  0.972 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        3.42800    0.04804  71.359  < 2e-16 ***
## Speciesversicolor -0.65800    0.06794  -9.685  < 2e-16 ***
## Speciesvirginica  -0.45400    0.06794  -6.683 4.54e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3397 on 147 degrees of freedom
## Multiple R-squared:  0.4008, Adjusted R-squared:  0.3926 
## F-statistic: 49.16 on 2 and 147 DF,  p-value: < 2.2e-16

Ensuite, on effectue le test post-hoc Tukey:

compSW <- aov(Sepal.Width ~ Species, data = iris)
TukeyHSD(compSW)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Sepal.Width ~ Species, data = iris)
## 
## $Species
##                        diff         lwr        upr     p adj
## versicolor-setosa    -0.658 -0.81885528 -0.4971447 0.0000000
## virginica-setosa     -0.454 -0.61485528 -0.2931447 0.0000000
## virginica-versicolor  0.204  0.04314472  0.3648553 0.0087802
summary(compSW)
##              Df Sum Sq Mean Sq F value Pr(>F)    
## Species       2  11.35   5.672   49.16 <2e-16 ***
## Residuals   147  16.96   0.115                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

On peut ensuite construire notre figure. D’abord, la fonction ggplot(). On spécifie le jeu de données, puis l’argument aes(), ici les axes x et y:

ggplot(iris, aes(x = Species, y = Sepal.Width))

Où sont les données? Il faut les ajouter!

C’est comme ça que ggplot() fonctionne. On crée la base de notre figure, puis on y ajoute les données à l’aide des fonctions geom_*. Par exemple, geom_points nous permet d’ajouter des données sous forme de scatterplots (donc des points), tandis que geom_line nous permet d’ajouter une ligne. Dans notre cas, nous voulons montrer nos données sous forme de boxplot, donc nous utilisons geom_boxplot. Il ne faut pas oublier d’ajouter un + après chaque fonction:

ggplot(iris, aes(x = Species, y = Sepal.Width)) +
geom_boxplot()

Pour ajouter/modifier d’autres éléments, il faut ajouter des couches. Ainsi, on peut changer les titres d’axes (labs pour labels), changer le thème de notre graphique (theme_bw est plus minimaliste), et avec la fonction theme(), retirer la grille et centrer le titre:

ggplot(iris, aes(x = Species, y = Sepal.Width)) +
geom_boxplot()+
labs(x = "Espèce", 
       y = "Largeur des sépales (cm)",
       title = "Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(), 
      plot.title = element_text(hjust = 0.5))

On peut aussi décider de changer les couleurs. Ce n’est pas toujours nécessaire ou pertinent, si ça n’ajoute pas d’information à la figure. Il faut faire attention de ne pas surcharger la figure avec le design esthétique, l’important c’est de focuser sur le message que l’on veut transmettre avec notre figure, et ne pas ajouter d’éléments qui sont distrayants. Mais pour l’exercice et se familiariser avec les paramètres graphiques, essayons:

ggplot(iris, aes(x = Species, y = Sepal.Width)) +
geom_boxplot(fill = "salmon", 
             color = "red")+
labs(x="Espèce", 
     y="Largeur des sépales (mm)",
     title="Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(), 
      plot.title = element_text(hjust = 0.5)) 

On a donc changé le fill (intérieur) et color (trait) de nos boxplot, à même la fonction boxplot(). On pourrait aussi changer la couleur des boîtes en fonction de l’espèce d’iris (qui est aussi notre variable x). Dans ce cas, il faut ajouter un argument aes() à la fonction boxplot(), pour assigner notre variable à un aesthetic. Ce sera notre variable Species:

#Fill:

ggplot(iris, aes(x = Species, y = Sepal.Width)) +
geom_boxplot(aes(fill = Species), 
             show.legend = FALSE)+
labs(x = "Espèce", 
     y = "Largeur des sépales (mm)",
     title = "Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(), 
      plot.title = element_text(hjust = 0.5)) 

#Color:
ggplot(iris, aes(x = Species, y = Sepal.Width)) +
geom_boxplot(aes(color = Species), 
             show.legend = FALSE)+
labs(x = "Espèce", 
     y = "Largeur des sépales (mm)",
     title = "Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(), 
      plot.title = element_text(hjust = 0.5)) 

Il y a aussi des palettes de couleurs prédéfinies dans R, qu’on peut utiliser avec la fonction scale_fill_brewer (j’ai choisi la palette “Accent”):

ggplot(iris, aes(x = Species, y = Sepal.Width)) +
geom_boxplot(aes(fill = Species), 
             show.legend = FALSE)+
scale_fill_brewer(palette = "Accent")+
labs(x = "Espèce", 
     y = "Largeur des sépales (mm)",
     title = "Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(), 
      plot.title = element_text(hjust = 0.5))

#Qu'est-ce qui arrive si on change Fill pour Color, dans boxplot()? Il faut aussi changer scale_fill_brewer() pour scale_color_brewer():

ggplot(iris, aes(x = Species, y=Sepal.Width)) +
geom_boxplot(aes(color = Species), 
             show.legend = FALSE)+
scale_color_brewer(palette = "Accent")+
labs(x = "Espèce", 
     y = "Largeur des sépales (mm)",
     title = "Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(), 
      plot.title = element_text(hjust = 0.5))

Puisque ggplot() fonctionne par couche, on peut aussi choisir d’ajouter les données brutes à notre figure. Par exemple, on ajoute avant geom_boxplot() une autre couche de données avec geom_point(), dans lequel on peut spécifier, encore une fois, l’argument aes() pour assigner la variable qui sera colorée. L’argument position définit le niveau de “jitter”, c’est-à-dire des points un peu “éparpillés” sur l’axe x. Ensuite, on pourrait ajouter l’argument alpha=0.5 dans la fonction boxplot(), pour les rendre semi-transparents:

ggplot(iris, aes(x = Species, y = Sepal.Width)) +
geom_point(aes(color = Species), 
           position = position_jitter(0.2), 
           shape = 16, 
           size = 2.5, 
           show.legend = FALSE)+
scale_color_brewer(palette = "Accent")+
geom_boxplot(alpha = 0.5)+
labs(x = "Espèce", 
     y = "Largeur des sépales (mm)", 
     title = "Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(), 
      plot.title = element_text(hjust = 0.5))

On peut aussi ajouter des lettres pour montrer nos résultats de comparaison (voir https://statdoe.com/one-way-anova-and-box-plot-in-r/ pour la méthode). En gros, la fonction multcompLetters4 assigne les lettres en fonction des groupes du test de Tukey, puis on constuit une table qui contient ces lettres (LettresSW dans l’exemple), la variable Species, ainsi que la position de la lettre (en haut du quartile pour chaque boîte/valeur de Species).

Deux nouvelles libraries sont nécessaires.

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(multcompView)

compSW <- aov(Sepal.Width ~ Species, data = iris)
tukeySW <- TukeyHSD(compSW)
print(tukeySW)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Sepal.Width ~ Species, data = iris)
## 
## $Species
##                        diff         lwr        upr     p adj
## versicolor-setosa    -0.658 -0.81885528 -0.4971447 0.0000000
## virginica-setosa     -0.454 -0.61485528 -0.2931447 0.0000000
## virginica-versicolor  0.204  0.04314472  0.3648553 0.0087802
cld <- multcompLetters4(compSW, tukeySW)

LettresSW <- group_by(iris, Species) %>%
  summarise(mean = mean(Sepal.Width), quant = quantile(Sepal.Width, probs = 0.75)) %>%
  arrange(desc(mean))

cld <- as.data.frame.list(cld$Species)
LettresSW$cld <- cld$Letters

print(LettresSW)
## # A tibble: 3 × 4
##   Species     mean quant cld  
##   <fct>      <dbl> <dbl> <chr>
## 1 setosa      3.43  3.68 a    
## 2 virginica   2.97  3.18 b    
## 3 versicolor  2.77  3    c

On ajoute nos lettres au boxplot avec la fonction geom_text(), avec comme argument aes() la table que l’on vient de créer. On peut aussi ajouter la statistique de test avec la fonction annotate(), en utilisant les coordonnées x et y dans notre figure.

ggplot(iris, aes(x = Species, y = Sepal.Width)) +
geom_point(aes(color = Species), 
           position = position_jitter(0.1), 
           shape = 16, 
           size = 2.5, 
           show.legend = FALSE)+
scale_color_brewer(palette = "Accent")+
geom_boxplot(alpha = 0.5)+
labs(x = "Espèce", 
     y = "Largeur des sépales (mm)",
     title = "Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(), 
      plot.title = element_text(hjust = 0.5))+
geom_text(data = LettresSW, 
          aes(x = Species, y = quant, label = cld), 
          size = 5, 
          vjust=-1, hjust =-5)+
annotate("text", x = 3, y = 4.5, label = "ANOVA, F(2)=49.16, p<0.001")

Pour bien faire, on pourrait même changer les étiquettes en x. Avec la fonction scale_x_discrete, j’ajoute “I.” pour iris à chaque étiquette et je les mets en italique, puisque ce sont des noms latins, avec l’argument axis.text.x dans la fonction theme().

ggplot(iris, aes(x = Species, y = Sepal.Width)) +
geom_point(aes(color = Species), 
           position = position_jitter(0.1), 
           shape = 16, 
           size = 2.5, show.legend = FALSE)+
scale_color_brewer(palette = "Accent")+
geom_boxplot(alpha = 0.5)+
labs(x = "Espèce", 
     y = "Largeur des sépales (mm)",
     title = "Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(), 
      plot.title = element_text(hjust = 0.5), 
      axis.text.x = element_text(face="italic"))+
geom_text(data = LettresSW, 
          aes(x = Species, y = quant, label = cld), 
          size = 5, 
          vjust=-1, hjust =-5)+
annotate("text", x = 3, y = 4.5, label = "ANOVA, F(2)=49.16, p<0.001")+
scale_x_discrete(labels = c("I. setosa", "I. versicolor", "I. virginica"))

Régression linéaire


GLM

Essayons avec un type d’analyse différent, les GLM binomiaux. Nous allons utiliser le jeu de données “survie”. Encore une fois, référez-vous à l’atelier sur les modèles GLM pour bien comprendre les analyses, j’y vais plus grossièrement ici.

Nous aurons besoin de:

library(ggeffects)

D’abord, télécharger le jeu de données “survie.csv”, puis l’ouvrir dans R (*le jeu provient du cours BIO8092).

surv <- read.delim2('survie.csv', header = TRUE, sep = ';', dec = '.')
head(surv)
##       masse      age sex survie
## 1 103.31373 60.43044   0      1
## 2 100.38447 63.11296   1      1
## 3 101.89734 60.85186   1      1
## 4  98.54894 62.10991   0      0
## 5 101.82275 59.57312   0      1
## 6 102.47745 56.40091   1      1

On a donc des données de masse, d’âge et de sexe (sous forme 0/1, mâle=1 et femelle=0), ainsi que ce qui sera notre variable réponse, la survie d’une bibitte quelconque. Celle-ci est aussi notée 0/1, car évidemment, c’est une variable binaire. Puisque celle-ci ne suit pas une loi normale, nous devons utiliser un GLM.

Prenez bien le temps de comprendre comment fonctionne la régression logistique. L’interprétation des coefficients n’est pas aussi directe ou intuitive que pour la régression linéaire. Comme vous vous souvenez sûrement, les paramètres sont dans une échelle différente que l’échelle de la “réalité”: la pente nous donne plutôt le changement dans le logit de la probabilité, pour un changement d’une unité de la variable explicative.

Rappel: Les GLM fonctionnent avec ce qu’on appelle une fonction de lien (“logit”, dans le cas d’une variable binaire - il en existe d’autre pour les autres distributions, comme celle de poisson ou gamma). Cette transformatiom permet de passer de l’échelle des valeurs observées (soit 0 et 1), à une échelle qui est linéaire, qui nous permet de “linéariser” un modèle à variable réponse binaire. Ainsi: Les valeurs observées de y (0 et 1) sont transformées avec la fonction de lien f() pour être sur la même échelle que le prédicteur linéaire (a+bxi). Les paramètres du modèle, nécessairement, sont dans cette échelle transformée. C’est pour cela qu’on ne peut pas calculer directement les valeurs prédites, comme avec la régression linéaire. Plutôt, les valeurs prédites par le modèle sont obtenues en appliquant l’inverse de la transformation, donc l’inverse de la fonction de lien (f−1), au prédicteur linéaire. Ainsi, Probabilité(succès)= logit−1(a+bx).

On peut convertir les paramètres du modèle logistique en ratios des cotes, en calculant e^-b, b étant la valeur d’un paramètre donné dans un modèle (voir https://statisticsbyjim.com/probability/odds-ratio/). Ce ratio des cotes nous indique le changement de cote associé à l’augmentation d’une unité de notre variable indépendante, par exemple, l’âge (s’il s’agit d’une autre variable binaire, par exemple le sexe, l’augmentation d’unité signifie qu’on passe de zéro, soit les femelles, à 1, soit les mâles − l’ordonnée à l’origine représente le 0 de toutes les variables).

En gros,

Si le ratio des cotes est >1, la probabilité de succès (dans notre cas, la survie) augmente avec l’augmentation de notre variable dépendante.

Si le ratio des cotes est <1, la probabilité de succès (dans notre cas, la survie) diminue avec l’augmentation de notre variable dépendante.

Calculer le ratio des cotes e^-b nous permet d’avoir une idée générale de l’effet d’une variable indépendante sur notre variable réponse.

Mais ce qui vous intéresse dans le cadre de cet atelier, c’est de visualiser les effets desdites variables!

Pour simplifier les calculs, on peut utiliser la fonction ggpredict() du package ggeffects. Celle-ci nous calcule les probabilités de survie (ou de succès) ainsi qu’un intervalle de confiance. On peut ensuite utiliser ces probabilités pour réaliser une figure avec ggplot().

Reprenons notre jeu de données.

head(surv)
##       masse      age sex survie
## 1 103.31373 60.43044   0      1
## 2 100.38447 63.11296   1      1
## 3 101.89734 60.85186   1      1
## 4  98.54894 62.10991   0      0
## 5 101.82275 59.57312   0      1
## 6 102.47745 56.40091   1      1

On peut commencer avec la question suivante: quel est l’effet de la masse sur la probabilité de survie?

Pour faire notre modèle, on utilise la fonction glm, qui ressemble à la fonction lm que nous avons utilisée plus tôt. Avec glm, on ajoute l’argument family = binomial, qui spécifie la distribution avec laquelle nous travaillons, soit la distribution binomiale.

modele.masse <- with(surv, glm(survie ~ masse, family = binomial))

summary(modele.masse)
## 
## Call:
## glm(formula = survie ~ masse, family = binomial)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.0758   0.2901   0.3887   0.5126   1.4173  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept) -45.9643    21.0574  -2.183   0.0290 *
## masse         0.4726     0.2095   2.256   0.0241 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 39.880  on 47  degrees of freedom
## Residual deviance: 34.075  on 46  degrees of freedom
## AIC: 38.075
## 
## Number of Fisher Scoring iterations: 5
#Quel serait le ratio des cotes?

b.masse<-modele.masse$coefficients['masse']

exp(b.masse)
##    masse 
## 1.604189

La masse aurait donc un effet positif sur la probabilité de survie. Essayons de visualiser ce résultat.

J’utilise la fonction ggpredict du package ggeffects. On spéficie d’abord le modèle à prédire, ensuite la variable que l’on souhaite visualiser en x.

predictions.masse<-ggpredict(modele.masse, "masse")
## Data were 'prettified'. Consider using `terms="masse [all]"` to get
##   smooth plots.
predictions.masse
## # Predicted probabilities of survie
## 
## masse | Predicted |       95% CI
## --------------------------------
##    96 |      0.36 | [0.07, 0.81]
##    97 |      0.47 | [0.14, 0.83]
##    99 |      0.70 | [0.44, 0.87]
##   100 |      0.79 | [0.60, 0.90]
##   101 |      0.85 | [0.71, 0.93]
##   102 |      0.90 | [0.76, 0.96]
##   103 |      0.94 | [0.80, 0.98]
##   106 |      0.98 | [0.85, 1.00]
## 
## Not all rows are shown in the ouput. Use `print(..., n = Inf)` to show
##   all rows.

Dans le récapitulatif, on voit l’étendue des valeurs de masse, ainsi que la probabilité de survie prédite pour chaque valeur de masse. De plus, un intervalle de confiance à 95% a été calculé. On a donc ce qu’il faut pour réaliser notre figure.

Nous utiliserons comme source de données, la table de prédictions que nous venons de sortir.

plot.masse <- ggplot()+
geom_line(data = predictions.masse, 
          aes(x = x, y = predicted))

plot.masse

On voit que la valeur maximale de l’axe y est de 1, qui correspondrait à une probabilité de 100%. Il est important d’ajouter notre invervalle de confiance, avec geom_ribbon. On utilise à nouveau notre table de prédictions, puisqu’elle contient les intervalles de confiance. Avec aes(ymin = conf.low, ymax = conf.high), on délimite les limites haut et bas de notre ruban. Ces limites correspondent aux valeurs d’intervalle calculées. Avec alpha = .2, le ruban est plus translucide:

plot.masse <- ggplot()+
geom_line(data = predictions.masse, 
          aes(x = x, y = predicted))+
geom_ribbon(data = predictions.masse, 
            show.legend = FALSE, 
            aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), 
            alpha = .2)

plot.masse

C’est bien, mais notre figure est pas mal tout nue!

On peut identifier nos axes + titre, ajouter le theme_bw(), enlever la grille + centrer le titre:

plot.masse <- ggplot()+
geom_line(data = predictions.masse, 
          aes(x = x, y = predicted))+
geom_ribbon(data = predictions.masse, 
            show.legend = FALSE, 
            aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), 
            alpha = .2)+
labs(x = "Masse (g)", 
     y = "Probabilité de survie",
     title = "Probabilité de survie de (?) en fonction de la masse")+
theme_bw()+
theme(panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(), 
      plot.title = element_text(hjust = 0.5))

plot.masse

Les valeurs de notre axe y sont de 0 à 1. On peut transformer notre axe pour qu’il soit exprimé en % sur 100, avec scale_y_continuous. De plus, ggplot() a choisi de commencer notre axe un peu en haut de zéro (même si c’est moins évident sur cette figure - retournez voir celle sans le ruban) car il optimise l’espace que prend nos données sur le graphique. Il serait plus juste de commencer notre axe à zéro. J’ajoute l’argument limits = c(0, 1) dans scale_y_continuous.

plot.masse <- ggplot()+
geom_line(data = predictions.masse, 
          aes(x = x, y = predicted))+
geom_ribbon(data = predictions.masse, 
            show.legend = FALSE, 
            aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), 
            alpha = .2)+
labs(x = "Masse (g)", 
     y = "Probabilité de survie", 
     title = "Probabilité de survie de (?) en fonction de la masse")+
theme_bw()+
theme(panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(), 
      plot.title = element_text(hjust = 0.5))+
scale_y_continuous(labels = scales::percent, 
                   limits = c(0, 1))

plot.masse

Si on voulait avoir moins de ticks sur l’axe des y (peut rendre la figure plus facilement lisible si elle est en petit format), on peut ajouter breaks = seq(0, 1, by = 0.5) dans scale_y_continuous, pour que les ticks soient à intervalle de 0.5:

plot.masse <- ggplot()+
geom_line(data = predictions.masse, 
          aes(x = x, y = predicted))+
geom_ribbon(data = predictions.masse, 
            show.legend = FALSE, 
            aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), 
            alpha = .2)+
labs(x = "Masse (g)", 
     y = "Probabilité de survie",
     title = "Probabilité de survie de (?) en fonction de la masse")+
theme_bw()+
theme(panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(), 
      plot.title = element_text(hjust = 0.5))+
scale_y_continuous(labels = scales::percent, 
                   limits = c(0, 1), 
                   breaks = seq(0, 1, by = 0.5))

plot.masse

Terminons la figure avec une jolie couleur, tiens: Dans geom_line(), j’ajoute color=group à l’argument aes(). On indique ainsi que le paramètre graphique de la couleur du trait sera associée aux “groupes”. Pour l’instant, on n’a qu’un seul groupe, donc ça peut sembler abstrait. Par contre, c’est ce qui nous permet, avec scale_color_manual(), de changer la couleur de notre ligne.

Même chose pour le ruban: dans geom_ribbon(), j’ajoute fill=group à l’argument aes(), pour assigner le remplissage du ruban aux “groupes”.

plot.masse <- ggplot()+
geom_line(data = predictions.masse, 
          aes(x = x, y = predicted, color = group), 
          show.legend = FALSE)+
scale_color_manual(values = c("#D28D32"))+
geom_ribbon(data = predictions.masse, 
            show.legend = FALSE, 
            aes(x = x, y = predicted, fill = group, ymin = conf.low, ymax = conf.high), 
            alpha = .2)+
scale_fill_manual(values = c("#D28D32"))+
labs(x = "Masse (g)", 
     y = "Probabilité de survie", 
     title = "Probabilité de survie de (?) en fonction de la masse")+
theme_bw()+
theme(panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(), 
      plot.title = element_text(hjust = 0.5))+
scale_y_continuous(labels = scales::percent, 
                   limits = c(0, 1), 
                   breaks = seq(0, 1, by = 0.5))

plot.masse

Remarquez que j’utilise color pour la ligne, et fill pour le ruban, comme nous avions abordé plus tôt avec les points dans l’ANOVA. En effet, on ne peut pas “fill” une ligne; en revanche, on peut “color” un ruban… Mais c’est sa bordure qui sera colorée:

plot.masse <- ggplot()+
geom_line(data = predictions.masse, 
          aes(x = x, y = predicted, color = group), 
          show.legend = FALSE)+
scale_color_manual(values = c("#D28D32"))+
geom_ribbon(data = predictions.masse, 
            show.legend = FALSE, 
            aes(x = x, y = predicted, color = group, ymin = conf.low, ymax = conf.high), 
            alpha = .2)+
theme_bw()+
theme(panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(), 
      plot.title = element_text(hjust = 0.5))

plot.masse

On pourrait aussi ajouter une couche de données contenant les données brutes, avec geom_point. Toutefois, on souhaite afficher les données brutes issues de surv, et non les résultats de notre modèle predictions.masse. Donc comme avec notre ANOVA, on assigne nos aes() dans geom_point(data=surv, aes(x=masse, y=survie)). J’utilise position = position_jitter(), car cela permet de mieux voir les points qui sont collés tous à la même valeur, soit 0 et 1 (aussi, détail, mais pour cette figure, je dois aussi ajouter height = 0.01 pour réduire le niveau de jitter appliqué en y, sinon le jitter s’énerve un peu trop et remplit tout l’espace - je crois que c’est lié au fait que notre variable n’a que deux valeurs…). De plus, vu que ça fait “déborder” les points un peu en dessous de 0 et un peu au dessus de 1, j’élargis les limites de l’axe y avec limits = c(-0.02, 1.02)) dans scale_y_continuous().

On termine en changeant la couleur des points et en augmentant la taille du trait avec linewidth=1.1 dans geom_line.

plot.masse <- ggplot()+
geom_line(data = predictions.masse, 
          aes(x = x, y = predicted, color = group), 
          show.legend = FALSE, 
          linewidth = 1.1)+
scale_color_manual(values = c("#D28D32"))+
geom_ribbon(data = predictions.masse, 
            show.legend = FALSE, 
            aes(x = x, y = predicted, fill = group, ymin = conf.low, ymax = conf.high), 
            alpha = .2)+
scale_fill_manual(values = c("#D28D32"))+
geom_point(data = surv, 
           aes(x = masse, y = survie), 
           alpha = .4, 
           color = "cornsilk4", 
           position = position_jitter(0.1, height = 0.01))+
labs(x = "Masse (g)", 
     y = "Probabilité de survie",
     title = "Probabilité de survie de (?) en fonction de la masse")+
theme_bw()+
theme(panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(), 
      plot.title = element_text(hjust = 0.5))+
scale_y_continuous(labels = scales::percent, 
                   limits = c(-0.02, 1.02))

plot.masse


Dans le jeu de données, on retrouve aussi la variable sexe:

head(surv)
##       masse      age sex survie
## 1 103.31373 60.43044   0      1
## 2 100.38447 63.11296   1      1
## 3 101.89734 60.85186   1      1
## 4  98.54894 62.10991   0      0
## 5 101.82275 59.57312   0      1
## 6 102.47745 56.40091   1      1

On pourrait donc faire un deuxième modèle qui inclut l’effet du sexe de la bibitte dans sa probabilité de survie.

modele.masse.sexe <- with(surv, glm(survie ~ masse + sex, family = binomial))
summary(modele.masse.sexe)
## 
## Call:
## glm(formula = survie ~ masse + sex, family = binomial)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3327   0.1828   0.3944   0.5588   1.3727  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept) -60.3269    26.0796  -2.313   0.0207 *
## masse         0.6094     0.2574   2.368   0.0179 *
## sex           1.3184     1.0860   1.214   0.2247  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 39.880  on 47  degrees of freedom
## Residual deviance: 32.401  on 45  degrees of freedom
## AIC: 38.401
## 
## Number of Fisher Scoring iterations: 6

Il ne semble pas avoir d’effet significatif du sexe sur la probabilité de survie.

Mais, pour l’exercice, faisons comme si l’effet était significatif et qu’on voulait l’illustrer dans une figure!

On calcule d’abord nos valeurs prédites de survie:

predictions.masse.sexe.test <- ggpredict(modele.masse.sexe, c("masse"))
## Data were 'prettified'. Consider using `terms="masse [all]"` to get
##   smooth plots.
predictions.masse.sexe.test
## # Predicted probabilities of survie
## 
## masse | Predicted |       95% CI
## --------------------------------
##    96 |      0.24 | [0.03, 0.76]
##    97 |      0.36 | [0.08, 0.78]
##    99 |      0.66 | [0.40, 0.85]
##   100 |      0.78 | [0.59, 0.90]
##   101 |      0.87 | [0.71, 0.95]
##   102 |      0.92 | [0.77, 0.98]
##   103 |      0.96 | [0.80, 0.99]
##   106 |      0.99 | [0.87, 1.00]
## 
## Adjusted for:
## * sex = 0.50
## 
## Not all rows are shown in the ouput. Use `print(..., n = Inf)` to show
##   all rows.

Remarquez dans le récapitulatif: Adjusted for: sex = 0.50

Rappelez-vous comment est construit un modèle GLM:

Les différents termes (paramètres β)*(la valeur de x) sont additionnés pour obtenir la variable réponse.

Donc, quand on lui demande de calculer la probabilité de survie en fonction de la masse, la fonction ggpredit() utilise par défaut la valeur moyenne de l’autre variable. Il set donc la valeur de sexe à 0,5. On obtient comme figure:

plot.masse.sexe.test <- ggplot()+
geom_line(data = predictions.masse.sexe.test, 
          aes(x = x, y = predicted, color = group), 
          show.legend = FALSE, 
          linewidth=1.1)+
labs(title = "Probabilité de survie ajustée pour sexe = 0.5")+
theme(plot.title = element_text(hjust = 0.5))

plot.masse.sexe.test

Dans notre cas, par contre, il serait plus intéressant de pouvour visualiser cette différence.

Pour cela, il faut inclure la variable sexe dans la fonction ggpredict:

predictions.masse.sexe <- ggpredict(modele.masse.sexe, c("masse","sex"))
## Data were 'prettified'. Consider using `terms="masse [all]"` to get
##   smooth plots.
predictions.masse.sexe
## # Predicted probabilities of survie
## 
## # sex = 0
## 
## masse | Predicted |       95% CI
## --------------------------------
##    96 |      0.14 | [0.01, 0.75]
##    98 |      0.35 | [0.07, 0.80]
##    99 |      0.50 | [0.17, 0.84]
##   101 |      0.77 | [0.52, 0.91]
##   103 |      0.92 | [0.73, 0.98]
##   106 |      0.99 | [0.84, 1.00]
## 
## # sex = 1
## 
## masse | Predicted |       95% CI
## --------------------------------
##    96 |      0.38 | [0.07, 0.83]
##    98 |      0.67 | [0.31, 0.90]
##    99 |      0.79 | [0.48, 0.94]
##   101 |      0.93 | [0.70, 0.99]
##   103 |      0.98 | [0.79, 1.00]
##   106 |      1.00 | [0.86, 1.00]
## 
## Not all rows are shown in the ouput. Use `print(..., n = Inf)` to show
##   all rows.

Cette fois-ci, nous obtenons deux tables de valeurs prédites - une pour sex = 0 (femelles) et une pour sex = 1 (mâles).

Dans le précédent modèle, nous avions défini l’aes(color=group), pour colorer notre ligne. La notion de groupe était plus abstraite, car il n’y avait en réalité qu’un seul groupe. Par contre, ici, dans plot.masse.sexe, nous en avons deux - les femelles et les mâles. Donc, quand on assigne color aux group, nous colorons les deux groupes de couleurs différentes:

plot.masse.sexe <- ggplot()+
geom_line(data = predictions.masse.sexe, 
          aes(x = x, y = predicted, color = group), 
          show.legend = TRUE, 
          linewidth=1.1)

plot.masse.sexe

On doit ajouter nos intervalles de confiance:

plot.masse.sexe <- ggplot()+
geom_line(data = predictions.masse.sexe, 
          aes(x = x, y = predicted, color = group), 
          show.legend = TRUE)+
geom_ribbon(data = predictions.masse.sexe, 
            show.legend = FALSE, 
            aes(x = x, y = predicted, fill = group, ymin = conf.low, ymax = conf.high), 
            alpha = 0.2)

plot.masse.sexe

La légende doit aussi être modifiée. scale_color_manual traite la couleur du trait formé par geom_line(), et scale_fill_manual, le remplissage du ruban geom_ribbon. En même temps, avec ces mêmes fonctions, on peut spécifier le titre de la légende (on n’écrit pas “Légende”!), le nom des groupes (la variable sex est codée 0 / 1, donc on doit ajouter manuellement ces noms). On peut aussi changer les couleurs.

plot.masse.sexe <- ggplot()+
geom_line(data = predictions.masse.sexe, 
          aes(x = x, y = predicted, color = group), 
          show.legend = TRUE)+
geom_ribbon(data = predictions.masse.sexe, 
            show.legend = TRUE, 
            aes(x = x, y = predicted, fill = group, ymin = conf.low, ymax = conf.high), 
            alpha = 0.2)+
scale_color_manual(name = "Sexe", 
                   labels = c("Femelle", "Mâles"), 
                   values = c("0" = "#D28D32", "1" = "darkblue"))+
scale_fill_manual(name = "Sexe", 
                  labels = c("Femelle", "Mâles"), 
                  values = c("0" = "#D28D32", "1" = "darkblue"))

plot.masse.sexe

Il est important que les noms et les aes() soient constants pour chaque ligne de notre figure ggplot(). Sinon, des légendes supplémentaires sont créées:

plot.masse.sexe <- ggplot()+
geom_line(data = predictions.masse.sexe, 
          aes(x = x, y = predicted, color = group), 
          show.legend = TRUE)+
geom_ribbon(data = predictions.masse.sexe, 
            show.legend = TRUE, 
            aes(x = x, y = predicted, fill = group, ymin = conf.low, ymax = conf.high), 
            alpha = 0.2)+
scale_color_manual(name = "Sexe", 
                   labels = c("Femelle", "Mâles"), 
                   values = c("0" = "#D28D32", "1" = "darkblue"))+
scale_fill_manual(name = "Coucou", 
                  labels = c("Femelle", "Mâles"), 
                  values = c("0" = "#D28D32", "1" = "darkblue"))

plot.masse.sexe

On peut ajouter les éléments graphiques que nous avions codés pour la figure de notre premier modèle: geom_line(linewidth=1.1) labs() theme_bw() theme() scale_y_continuous(labels = scales::percent, limits = c(-0.02, 1.02))

plot.masse.sexe <- ggplot()+
geom_line(data = predictions.masse.sexe, aes(x=x, y = predicted, color = group), 
          show.legend = TRUE, 
          linewidth = 1.1)+
geom_ribbon(data = predictions.masse.sexe, 
            show.legend = TRUE, 
            aes(x = x, y = predicted, fill = group, ymin = conf.low, ymax = conf.high), 
            alpha = 0.2)+
scale_color_manual(name = "Sexe", 
                   labels = c("Femelle", "Mâles"), 
                   values = c("0" = "#D28D32", "1" = "darkblue"))+
scale_fill_manual(name = "Sexe", 
                  labels = c("Femelle", "Mâles"), 
                  values = c("0" = "#D28D32", "1" = "darkblue"))+
labs(x = "Masse (g)", 
     y = "Probabilité de survie",
     title = "Probabilité de survie de (?) en fonction de la masse, pour les individus mâles et femelles")+
theme_bw()+
theme(panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(), 
      plot.title = element_text(hjust = 0.5))+
scale_y_continuous(labels = scales::percent, 
                   limits = c(-0.02, 1.02))

plot.masse.sexe

Le titre est devenu plus long et pourrait déborder de la figure (cela dépendra des paramètres choisis au moment d’exporter la figure). On peut utiliser \n pour indiquer précisément l’endroit où couper notre texte:

plot.masse.sexe <- ggplot()+
geom_line(data = predictions.masse.sexe, 
          aes(x = x, y = predicted, color = group), 
          show.legend = TRUE, 
          linewidth = 1.1)+
geom_ribbon(data = predictions.masse.sexe, 
            show.legend = TRUE, 
            aes(x = x, y = predicted, fill = group, ymin = conf.low, ymax = conf.high), 
            alpha = 0.2)+
scale_color_manual(name = "Sexe", 
                   labels = c("Femelle", "Mâles"), 
                   values = c("0" = "#D28D32", "1" = "darkblue"))+
scale_fill_manual(name = "Sexe", 
                  labels = c("Femelle", "Mâles"), 
                  values = c("0" = "#D28D32", "1" = "darkblue"))+
labs(x = "Masse (g)", 
     y = "Probabilité de survie",
     title = "Probabilité de survie de (?) en fonction de la\nmasse, pour les individus mâles et femelles")+
theme_bw()+
theme(panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(), 
      plot.title = element_text(hjust = 0.5))+
scale_y_continuous(labels = scales::percent, 
                   limits = c(-0.02, 1.02))

plot.masse.sexe

On peut ajouter nos données brutes pour terminer la figure.

Je dois faire une dernière étape pour afficher correctement les données brutes.

surv$sex<-as.factor(surv$sex)

Sinon, il traite sex comme étant une variable continue (ce qu’elle est), et on ne peut pas appliquer les paramètres graphiques correctement!

(Idéalement, faire cela avant l’analyse…!)

plot.masse.sexe <- ggplot()+
geom_line(data = predictions.masse.sexe, 
          aes(x = x, y = predicted, color = group), 
          show.legend = TRUE, 
          linewidth = 1.1)+
scale_color_manual(name = "Sexe", 
                   labels = c("Femelle", "Mâles"), 
                   values = c("0" = "#D28D32", "1" = "darkblue"))+
geom_ribbon(data = predictions.masse.sexe, 
            show.legend = TRUE, 
            aes(x = x, y = predicted, fill = group, ymin = conf.low, ymax = conf.high), 
            alpha = 0.2)+
scale_fill_manual(name = "Sexe", 
                  labels = c("Femelle", "Mâles"), 
                  values = c("0" = "#D28D32", "1" = "darkblue"))+
geom_point(data = surv, aes(x = masse, y = survie, col = sex), 
           alpha = .4, 
           position = position_jitter(0.1, height = 0.01))+
labs(x = "Masse (g)", 
     y = "Probabilité de survie",
     title = "Probabilité de survie de (?) en fonction de la\nmasse, pour les individus mâles et femelles")+
theme_bw()+
theme(panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(), 
      plot.title = element_text(hjust = 0.5))+
scale_y_continuous(labels = scales::percent, 
                   limits = c(-0.02, 1.02))+
scale_x_continuous(breaks = seq(96, 106, by = 2))

plot.masse.sexe


Si on avait une troisième variable dans notre modèle, par exemple l’âge:

head(surv)
##       masse      age sex survie
## 1 103.31373 60.43044   0      1
## 2 100.38447 63.11296   1      1
## 3 101.89734 60.85186   1      1
## 4  98.54894 62.10991   0      0
## 5 101.82275 59.57312   0      1
## 6 102.47745 56.40091   1      1

On peut l’ajouter à notre modèle:

modele.masse.sexe.age <- with(surv, glm(survie ~ masse + sex + age, family = binomial))

summary(modele.masse.sexe.age)
## 
## Call:
## glm(formula = survie ~ masse + sex + age, family = binomial)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4538   0.2002   0.3538   0.5734   1.2288  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept) -69.8689    30.8379  -2.266   0.0235 *
## masse         0.6183     0.2599   2.379   0.0173 *
## sex1          1.4678     1.1074   1.325   0.1850  
## age           0.1434     0.2213   0.648   0.5171  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 39.880  on 47  degrees of freedom
## Residual deviance: 31.972  on 44  degrees of freedom
## AIC: 39.972
## 
## Number of Fisher Scoring iterations: 6

Ok, cette variable n’est pas significative non plus. Mais encore une fois, essayons de la visualiser pour l’exercice.

On débute avec le calcul des prédictions, en ajoutant la variable age:

predictions.masse.sexe.age <- ggpredict(modele.masse.sexe.age, c("masse","sex","age"))
## Data were 'prettified'. Consider using `terms="masse [all]"` to get
##   smooth plots.
predictions.masse.sexe.age
## # Predicted probabilities of survie
## 
## # sex = 0
## # age = 57.64
## 
## masse | Predicted |       95% CI
## --------------------------------
##    96 |      0.10 | [0.00, 0.72]
##    99 |      0.40 | [0.08, 0.83]
##   101 |      0.70 | [0.33, 0.92]
##   106 |      0.98 | [0.76, 1.00]
## 
## # sex = 1
## # age = 57.64
## 
## masse | Predicted |       95% CI
## --------------------------------
##    96 |      0.32 | [0.05, 0.82]
##    99 |      0.75 | [0.38, 0.93]
##   101 |      0.91 | [0.63, 0.98]
##   106 |      1.00 | [0.84, 1.00]
## 
## # sex = 0
## # age = 59.98
## 
## masse | Predicted |       95% CI
## --------------------------------
##    96 |      0.13 | [0.01, 0.73]
##    99 |      0.49 | [0.16, 0.83]
##   101 |      0.77 | [0.51, 0.91]
##   106 |      0.99 | [0.83, 1.00]
## 
## # sex = 1
## # age = 59.98
## 
## masse | Predicted |       95% CI
## --------------------------------
##    96 |      0.39 | [0.07, 0.84]
##    99 |      0.80 | [0.49, 0.95]
##   101 |      0.93 | [0.71, 0.99]
##   106 |      1.00 | [0.87, 1.00]
## 
## # sex = 0
## # age = 62.32
## 
## masse | Predicted |       95% CI
## --------------------------------
##    96 |      0.17 | [0.01, 0.80]
##    99 |      0.57 | [0.18, 0.89]
##   101 |      0.82 | [0.51, 0.95]
##   106 |      0.99 | [0.84, 1.00]
## 
## # sex = 1
## # age = 62.32
## 
## masse | Predicted |       95% CI
## --------------------------------
##    96 |      0.47 | [0.07, 0.91]
##    99 |      0.85 | [0.45, 0.98]
##   101 |      0.95 | [0.68, 0.99]
##   106 |      1.00 | [0.87, 1.00]
## 
## Not all rows are shown in the ouput. Use `print(..., n = Inf)` to show
##   all rows.

Nous obtenons 6 tables de prédictions. Puisque age est une variable continue, ggpredict nous a choisi 3 valeurs de age représentatives de la distribution et a calculé un y prédit pour chacunes de ces valeurs et pour les 2 niveaux de sex.

Ici nous avons deux options. Puisque nous avons une troisième variable indépendante à représenter, on peut lui assigner un aes(), par exemple linetype.

Nous avons donc, dans nos tables de prédiction:

masse = x = x

survie = y = predicted

sex = color = group

Dans la table de prédiction, bien que ce ne soit pas affiché, la troisième variable indépendante est appelée facet, comme on peut remarquer dans

str(predictions.masse.sexe.age)
## Classes 'ggeffects' and 'data.frame':    66 obs. of  7 variables:
##  $ x        : int  96 96 96 96 96 96 97 97 97 97 ...
##  $ predicted: num  0.0958 0.1291 0.1717 0.315 0.3914 ...
##  $ std.error: num  1.63 1.48 1.51 1.16 1.08 ...
##  $ conf.low : num  0.00433 0.00803 0.01063 0.04535 0.07209 ...
##  $ conf.high: num  0.721 0.731 0.8 0.817 0.842 ...
##  $ group    : Factor w/ 2 levels "0","1": 1 1 1 2 2 2 1 1 1 2 ...
##  $ facet    : Factor w/ 3 levels "57.64","59.98",..: 1 2 3 1 2 3 1 2 3 1 ...
##  - attr(*, "numeric.facet")= logi TRUE
##  - attr(*, "legend.labels")= chr [1:2] "0" "1"
##  - attr(*, "x.is.factor")= chr "0"
##  - attr(*, "continuous.group")= logi FALSE
##  - attr(*, "rawdata")='data.frame':  48 obs. of  5 variables:
##   ..$ response: int [1:48] 1 1 1 0 1 1 1 1 1 1 ...
##   ..$ x       : num [1:48] 103.3 100.4 101.9 98.5 101.8 ...
##   ..$ group   : Factor w/ 2 levels "0","1": 1 2 2 1 1 2 1 2 1 2 ...
##   ..$ facet   : num [1:48] 60.4 63.1 60.9 62.1 59.6 ...
##   ..$ rowname : chr [1:48] "1" "2" "3" "4" ...
##  - attr(*, "title")= chr "Predicted probabilities of survie"
##  - attr(*, "x.title")= chr "masse"
##  - attr(*, "y.title")= chr "survie"
##  - attr(*, "legend.title")= chr "sex"
##  - attr(*, "constant.values")= Named list()
##  - attr(*, "terms")= chr [1:3] "masse" "sex" "age"
##  - attr(*, "original.terms")= chr [1:3] "masse" "sex" "age"
##  - attr(*, "at.list")=List of 3
##   ..$ masse: int [1:11] 96 97 98 99 100 101 102 103 104 105 ...
##   ..$ sex  : chr [1:2] "0" "1"
##   ..$ age  : num [1:3] 57.6 60 62.3
##  - attr(*, "ci.lvl")= num 0.95
##  - attr(*, "type")= chr "fe"
##  - attr(*, "response.name")= chr "survie"
##  - attr(*, "back.transform")= logi TRUE
##  - attr(*, "response.transform")= chr "survie"
##  - attr(*, "untransformed.predictions")= num [1:66] 0.0958 0.1291 0.1717 0.315 0.3914 ...
##  - attr(*, "family")= chr "binomial"
##  - attr(*, "link")= chr "logit"
##  - attr(*, "logistic")= chr "1"
##  - attr(*, "link_inverse")=function (eta)  
##  - attr(*, "link_function")=function (mu)  
##  - attr(*, "is.trial")= chr "0"
##  - attr(*, "fitfun")= chr "glm"
##  - attr(*, "verbose")= logi TRUE
##  - attr(*, "model.name")= chr "modele.masse.sexe.age"

Dans ce cas, on ajoute l’aes() linetype à geom_line assigné à facet:

plot.masse.sexe.age1 <- ggplot()+
geom_line(data = predictions.masse.sexe.age, 
          aes(x = x, y = predicted, color = group, linetype = facet), 
          show.legend = TRUE, 
          linewidth=1.1)

plot.masse.sexe.age1

On peut voir dans la légende les 3 valeurs de age qui ont été utilisées pour le calcul.

La deuxième option serait d’utiliser des facettes. Ceci crée des graphiques côte à côte avec différents niveaux d’une variable choisie. Dans ce cas, j’enlève l’aes() qu’on vient d’ajouter, et j’ajoute la fonction facet_wrap(~facet), qui crée 3 facettes avec notre “variable” facet:

plot.masse.sexe.age2 <- ggplot()+
geom_line(data = predictions.masse.sexe.age, 
          aes(x = x, y = predicted, color = group), 
          show.legend = TRUE, 
          linewidth = 1.1)+
facet_wrap(~ facet)

plot.masse.sexe.age2

On peut terminer ce graphique, avec le même code que les précédents. J’ajoute labeller dans facet_wrap pour changer le titre des étiquettes. J’ai arrondi l’âge pour les étiquettes.

plot.masse.sexe.age2 <- ggplot()+
geom_line(data = predictions.masse.sexe.age, 
          aes(x = x, y = predicted, color = group), 
          show.legend = TRUE, 
          linewidth = 1.1)+
scale_color_manual(name = "Sexe", 
                   labels = c("Femelle", "Mâles"), 
                   values = c("0" = "#D28D32", "1" = "darkblue"))+
geom_ribbon(data = predictions.masse.sexe, 
            show.legend = TRUE, 
            aes(x = x, y = predicted, fill = group, ymin = conf.low, ymax = conf.high), 
            alpha = 0.2)+
scale_fill_manual(name = "Sexe", 
                  labels = c("Femelle", "Mâles"), 
                  values = c("0" = "#D28D32", "1" = "darkblue"))+
geom_point(data = surv, 
           aes(x = masse, y = survie, col = sex), 
           alpha = .4, 
           position = position_jitter(0.1, height = 0.01))+
labs(x = "Masse (g)", 
     y = "Probabilité de survie", 
     title = "Probabilité de survie de (?) en fonction de la\nmasse et de l'âge, pour les individus mâles et femelles")+
theme_bw()+
theme(panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(), 
      plot.title = element_text(hjust = 0.5))+
scale_y_continuous(labels = scales::percent, 
                   limits = c(-0.02, 1.02))+
scale_x_continuous(breaks = seq(96, 106, by = 2))+
facet_wrap(~ facet, labeller = labeller(facet = c( "57.64" = "58 jours", "59.98" = "60 jours",  "62.32" = "62 jours"))) 

plot.masse.sexe.age2

Vous pouvez vous pratiquer en complétant la première option!

plot.masse.sexe.age1


Régression linéaire simple

Maintenant qu’on a fait un exemple plus compliqué, essayons d’illustrer un modèle linéaire simple.

Reprenons les données iris()

head(iris)
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa

On peut créer un modèle linéaire pour prédire Sepal.Length en fonction de Petal.Length:

Encore une fois, ce n’est qu’un exemple pour s’exercer avec ggplot(). Il ne faut pas oublier les étapes préliminaires, que je saute ici, avant de se lancer dans un modèle.

regr.SL <- lm(Sepal.Length ~ Petal.Length, data = iris)
summary(regr.SL)
## 
## Call:
## lm(formula = Sepal.Length ~ Petal.Length, data = iris)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.24675 -0.29657 -0.01515  0.27676  1.00269 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   4.30660    0.07839   54.94   <2e-16 ***
## Petal.Length  0.40892    0.01889   21.65   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4071 on 148 degrees of freedom
## Multiple R-squared:   0.76,  Adjusted R-squared:  0.7583 
## F-statistic: 468.6 on 1 and 148 DF,  p-value: < 2.2e-16

Enfin quelque chose de significatif! Super, nous pouvons représenter cettre relation.

Avec une régression linéaire simple, évidemment, nous n’avons pas besoin d’une fonction comme ggpredict pour calculer y.

On peut commencer par spécifier les aes() du jeu de données brutes directement dans ggplot(); ensuite, nous visualisons ces données avec geom_point; puis, geom_smooth(method = lm) (calcule, puis) ajoute la courbe de la régression.

plot.regr.SL <- ggplot(data = iris, aes(x = Petal.Length, y = Sepal.Length))+
geom_point()+
geom_smooth(method = lm, 
            se = TRUE)

plot.regr.SL
## `geom_smooth()` using formula = 'y ~ x'

On peut embellir notre graphique:

plot.regr.SL <- ggplot(data = iris, aes(x = Petal.Length, y = Sepal.Length))+
geom_point(col = "#8F1E5F")+
geom_smooth(method = lm, 
            se = TRUE, 
            col = "#2B091D", 
            fill = "#591130")+
labs(x = "Longueur du pétale (cm)", 
     y = "Longueur du sépale (mm)",
     title = "Longueur du pétale en fonction de la longueur du sépale chez les iris")+
theme_bw()+
theme(panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(), 
      plot.title = element_text(hjust = 0.5))

plot.regr.SL
## `geom_smooth()` using formula = 'y ~ x'

On retourne voir le résumé de notre modèle pour ajouter la statistique:

summary(regr.SL)
## 
## Call:
## lm(formula = Sepal.Length ~ Petal.Length, data = iris)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.24675 -0.29657 -0.01515  0.27676  1.00269 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   4.30660    0.07839   54.94   <2e-16 ***
## Petal.Length  0.40892    0.01889   21.65   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4071 on 148 degrees of freedom
## Multiple R-squared:   0.76,  Adjusted R-squared:  0.7583 
## F-statistic: 468.6 on 1 and 148 DF,  p-value: < 2.2e-16

Et on l’ajoute au graphique:

plot.regr.SL <- ggplot(data = iris, aes(x = Petal.Length, y = Sepal.Length))+
geom_point(col="#8F1E5F")+
geom_smooth(method = lm, 
            se = TRUE, col = "#2B091D", fill = "#591130")+
labs(x = "Longueur du pétale (cm)", 
     y = "Longueur du sépale (mm)",
     title = "Longueur du pétale en fonction de la longueur du sépale chez les iris")+
theme_bw()+
theme(panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(), 
      plot.title = element_text(hjust = 0.5))+
annotate("text", 
         x = 2.5, y = 7.75, 
         label = "R2 = 0.76, F(1, 148) = 468.6, p = < 0.001")

plot.regr.SL
## `geom_smooth()` using formula = 'y ~ x'

C’est important d’ajouter un intervalle de confiance, ainsi qu’un intervalle de prédiction. L’intervalle de confiance s’affiche avec geom_smooth(se = TRUE). Pour ce qui est de l’intervalle de prédiction, la fonction predict() calcule ces valeurs, puis colle ce nouveau dataframe à iris pour créer iris.int.pred. Nous utiliserons maintenant cette nouvelle table pour notre régression.

int.pred <- predict(regr.SL, interval = "prediction")
## Warning in predict.lm(regr.SL, interval = "prediction"): predictions on current data refer to _future_ responses
iris.int.pred <- cbind(iris, int.pred)
head(iris.int.pred)
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species      fit      lwr
## 1          5.1         3.5          1.4         0.2  setosa 4.879095 4.067202
## 2          4.9         3.0          1.4         0.2  setosa 4.879095 4.067202
## 3          4.7         3.2          1.3         0.2  setosa 4.838202 4.025897
## 4          4.6         3.1          1.5         0.2  setosa 4.919987 4.108491
## 5          5.0         3.6          1.4         0.2  setosa 4.879095 4.067202
## 6          5.4         3.9          1.7         0.4  setosa 5.001771 4.191017
##        upr
## 1 5.690987
## 2 5.690987
## 3 5.650508
## 4 5.731483
## 5 5.690987
## 6 5.812526

Dans geom_ribbon(aes(ymin = lwr, ymax = upr), les aes() permettent d’assigner les bornes supérieures et inférieures (en y) du ruban aux variables lwr et upr que nous venons de calculer et ajouter à notre table.

plot.regr.SL <- ggplot(data = iris.int.pred, aes(x = Petal.Length, y = Sepal.Length))+
geom_point(col = "#8F1E5F")+
geom_ribbon(aes(ymin = lwr, ymax = upr), 
            alpha = .2, 
            fill = "#591130")+
geom_smooth(method = lm, 
            se = TRUE, 
            col = "#2B091D", 
            fill = "#591130")+
labs(x = "Longueur du pétale (cm)", 
     y = "Longueur du sépale (mm)",
     title = "Longueur du pétale en fonction de la longueur du sépale chez les iris")+
theme_bw()+
theme(panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(), 
      plot.title = element_text(hjust = 0.5))+
annotate("text", 
         x = 2.5, y = 7.75, 
         label = "R2 = 0.76, F(1, 148) = 468.6, p = < 0.001")

plot.regr.SL
## `geom_smooth()` using formula = 'y ~ x'

Wow!


Cheat sheet ggplot()

Maintenant que vous êtes plus familiers avec ggplot(), voici un lien pour une cheat sheet pour ggplot():

https://rstudio.github.io/cheatsheets/html/data-visualization.html

Elle contient toute l’information sur les différents aes(), geom_*(), scale_*, ainsi que sur les thèmes, les ajustements de position, et les légendes, par exemple.


À considérer lors de la création de figures

Comme mentionné plus tôt, il faut être judicieux dans le choix de notre symbologie. En effet, les couleurs et les symboles ont souvent un sens qui nous est intuitif, et représentent ainsi de l’information - quelles couleurs choisiriez-vous pour illustrer une température chaude? Froide?

Ces choix de couleur contribuent à donner un sens à vos figures. Il existe des palettes de couleurs conçue pour illustrer des données, par exemple, d’élévation, ou de température. Certaines sont conçues pour maximiser le contraste entre les extrêmes des valeurs d’une variable. Il existe même une palette, “Okabe-Ito”, contenant des couleurs qui sont plus accessibles pour les gens ayant une variété de problèmes visuels. D’autres palettes sont simplement esthétiques.

On doit aussi penser au type de variable que nous souhaitons illustrer. Par exemple, il serait moins approprié d’utiliser la valeur (donc de pâle à foncé) d’une même couleur pour colorer des variables catégoriques qui ne sont pas ordonnées (pensez aux boxplots des espèces d’iris), puisque celle-ci traduit la notion d’ordre et de progression. Une palette avec des couleurs aléatoires comme celle que nous avons utilisée est mieux dans ce cas. La valeur serait toutefois appropriée pour des données de concentration, par exemple.

Toutefois, certaines publications scientifiques n’acceptent que les figures en noir et blanc; dans ce cas, on peut utiliser des symboles, des textures (par exemple, le grain peut être utilisé au lieu de la valeur), ou différentes lignes pointillées, par exemple.

De plus, assurez-vous que vos axes sont représentatifs de vos données.

Pour terminer, utilisez une fonction pour exporter vos figures. “Copier-coller” ne suffit pas, car la résolution ne sera pas assez bonne.

Voici un exemple:

ggsave(path = "votre_chemin_comme_dans_setwd()", "titre_de_votre_figures.jpeg", width = 6, height = 6, dpi=700)


Ordinations

Comme pour les analyses univariées, cet atelier est seulement axé sur la création de graphiques. Pour une introduction aux analyses multivariées (ordinations, PERMANOVA, Procrustes), consulter l’atelier 6.

Importer les données multivariées

library(vegan)
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.6-4
data(dune) # matrice de communautés
data(dune.env) # données environnementales

Ce jeu de données contient une matrice d’abondance (dune) d’espèces végétales avec leur classe de couverture pour 20 sites. Le dataframe (dune.env) contient 5 variables denvironnementales.

PCA

Préparation des données

Transformation Hellinger

dune.hel <- decostand(dune, method="hellinger")

Faire une PCA (avec la fonction prcomp cette fois-ci, car l’objet qu’on va créer peut être “lu” plus tard par la fonction ggord)

pca.dune <- prcomp(dune.hel)
summary(pca.dune)
## Importance of components:
##                           PC1    PC2    PC3    PC4     PC5     PC6     PC7
## Standard deviation     0.4122 0.3333 0.2374 0.2097 0.20242 0.17376 0.14506
## Proportion of Variance 0.3057 0.1998 0.1013 0.0791 0.07371 0.05431 0.03785
## Cumulative Proportion  0.3057 0.5054 0.6068 0.6859 0.75960 0.81391 0.85177
##                            PC8     PC9    PC10    PC11    PC12    PC13    PC14
## Standard deviation     0.13392 0.11745 0.10441 0.09995 0.09550 0.07985 0.07188
## Proportion of Variance 0.03226 0.02481 0.01961 0.01797 0.01641 0.01147 0.00929
## Cumulative Proportion  0.88403 0.90885 0.92846 0.94643 0.96283 0.97430 0.98360
##                           PC15    PC16    PC17    PC18    PC19      PC20
## Standard deviation     0.05424 0.05240 0.04063 0.03373 0.02536 1.816e-17
## Proportion of Variance 0.00529 0.00494 0.00297 0.00205 0.00116 0.000e+00
## Cumulative Proportion  0.98889 0.99383 0.99680 0.99884 1.00000 1.000e+00

Dans l’atelier 6, nous avons fait nos PCA avec la fonction rda. Lorsqu’on compare le sommaire des deux analyses, on voit que les deux fonctions donnent les mêmes résultats.

pca2.dune <- rda(dune.hel)
summary(pca2.dune)
## 
## Call:
## rda(X = dune.hel) 
## 
## Partitioning of variance:
##               Inertia Proportion
## Total          0.5559          1
## Unconstrained  0.5559          1
## 
## Eigenvalues, and their contribution to the variance 
## 
## Importance of components:
##                          PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Eigenvalue            0.1699 0.1111 0.05634 0.04397 0.04097 0.03019 0.02104
## Proportion Explained  0.3057 0.1998 0.10135 0.07910 0.07371 0.05431 0.03785
## Cumulative Proportion 0.3057 0.5054 0.60679 0.68590 0.75960 0.81391 0.85177
##                           PC8     PC9    PC10     PC11    PC12     PC13
## Eigenvalue            0.01794 0.01379 0.01090 0.009991 0.00912 0.006376
## Proportion Explained  0.03226 0.02481 0.01961 0.017972 0.01641 0.011469
## Cumulative Proportion 0.88403 0.90885 0.92846 0.946428 0.96283 0.974303
##                           PC14     PC15     PC16     PC17     PC18      PC19
## Eigenvalue            0.005166 0.002942 0.002745 0.001651 0.001138 0.0006432
## Proportion Explained  0.009293 0.005292 0.004939 0.002969 0.002047 0.0011570
## Cumulative Proportion 0.983597 0.988888 0.993827 0.996796 0.998843 1.0000000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores:  1.80276 
## 
## 
## Species scores
## 
##                 PC1       PC2       PC3       PC4        PC5        PC6
## Achimill -0.2241355  0.057048 -0.068674 -0.173281 -0.0498964 -0.0006022
## Agrostol  0.4142446 -0.180068  0.020078 -0.020756  0.0029054  0.0193453
## Airaprae -0.0350694  0.154644  0.116317 -0.039999 -0.1384374 -0.0376930
## Alopgeni  0.1386550 -0.318444  0.186419  0.001329  0.0030171 -0.0215368
## Anthodor -0.1962910  0.238579  0.095768 -0.179434 -0.0380470 -0.0471885
## Bellpere -0.1148332 -0.054380 -0.057640  0.064957  0.0010453  0.0896219
## Bromhord -0.1410948 -0.053134 -0.049586 -0.057365  0.0210053  0.0938974
## Chenalbu  0.0121539 -0.025919  0.036102 -0.015353 -0.0038179  0.0173353
## Cirsarve -0.0026599 -0.031975  0.004872  0.021813 -0.0182940  0.0064877
## Comapalu  0.1100812  0.048635 -0.061981 -0.028257  0.0102828  0.1064664
## Eleopalu  0.3796078  0.069099 -0.185862 -0.064110  0.0144632  0.0035959
## Elymrepe -0.1456860 -0.248314 -0.094249  0.010195 -0.1482516 -0.0563573
## Empenigr  0.0004895  0.057213  0.063829  0.026314 -0.0393467 -0.0066958
## Hyporadi -0.0559119  0.197829  0.139068  0.036814 -0.1427093 -0.0326104
## Juncarti  0.2478729 -0.012497 -0.092252  0.005525  0.0265396 -0.1569375
## Juncbufo  0.0190348 -0.121457  0.179346 -0.051491  0.0931930 -0.0211675
## Lolipere -0.3366259 -0.173707 -0.206664  0.154192  0.0092155 -0.0325721
## Planlanc -0.2547353  0.197554 -0.034482 -0.039980  0.1469090 -0.0256698
## Poaprat  -0.2808972 -0.142814 -0.082463  0.084918 -0.0502270 -0.0278725
## Poatriv  -0.1317430 -0.363370  0.059442 -0.137436  0.0581366 -0.0136469
## Ranuflam  0.2783344  0.024399 -0.081013 -0.052662  0.0009058  0.0274683
## Rumeacet -0.1079283 -0.024094  0.054027 -0.122507  0.2108960 -0.0801798
## Sagiproc  0.0409392 -0.082141  0.246528  0.127747 -0.0148906 -0.0086761
## Salirepe  0.0651846  0.156446  0.007210  0.136628 -0.0166452 -0.0581183
## Scorautu -0.0562239  0.158686  0.094988  0.095244  0.0265794  0.1029567
## Trifprat -0.1001443  0.028383 -0.024352 -0.095019  0.1389927 -0.0417313
## Trifrepe -0.0471194 -0.006999  0.052812  0.045812  0.1395303  0.2537096
## Vicilath -0.0541878  0.052774 -0.023104  0.114929  0.0381791  0.0343033
## Bracruta  0.0852752  0.107955 -0.003783  0.181534  0.2198174 -0.1398609
## Callcusp  0.2045115  0.057120 -0.104072 -0.065216 -0.0146101  0.0579472
## 
## 
## Site scores (weighted sums of species scores)
## 
##          PC1      PC2        PC3      PC4       PC5       PC6
## 1  -0.423423 -0.42754 -6.533e-01  0.01583 -0.761705 -0.463469
## 2  -0.356333 -0.36816 -1.991e-01 -0.08370 -0.355682  0.630146
## 3  -0.080900 -0.55215 -5.753e-02  0.27275 -0.217514  0.016724
## 4  -0.041005 -0.49292  7.510e-02  0.33626 -0.282018  0.100013
## 5  -0.458003  0.05220 -1.432e-01 -0.54110  0.351891 -0.117334
## 6  -0.389342  0.20337 -1.022e-01 -0.30549  0.758789 -0.245696
## 7  -0.451813  0.06864 -6.835e-02 -0.41821  0.585536 -0.138732
## 8   0.295113 -0.29375 -6.362e-02  0.19828  0.010903 -0.185034
## 9   0.035973 -0.49264  2.714e-01  0.07945  0.086166 -0.419257
## 10 -0.453094  0.13639 -2.391e-01 -0.12201  0.144355  0.462946
## 11 -0.273214  0.29633  3.931e-05  0.87663  0.126762  0.097096
## 12  0.246622 -0.33173  9.205e-01 -0.03527  0.493526 -0.017653
## 13  0.226908 -0.48390  6.740e-01 -0.28663 -0.071279  0.323642
## 14  0.570818  0.25288 -3.165e-01 -0.32869 -0.085439  1.196310
## 15  0.654415  0.28845 -3.733e-01  0.01035  0.196968  0.002255
## 16  0.720030 -0.10689 -2.727e-01 -0.35992  0.007984 -0.504518
## 17 -0.317465  0.75272  3.395e-01 -0.64284 -0.803236 -0.262492
## 18 -0.201131  0.39818 -2.007e-01  0.89875  0.365682  0.086314
## 19  0.006263  0.73205  8.167e-01  0.33669 -0.503443 -0.085674
## 20  0.689579  0.36848 -4.077e-01  0.09888 -0.048247 -0.475587

PERMANOVA. On va tester si la composition végétale diffère entre différents type d’aménagement (facteur; variable quantitative) des sites.

dune.dist <- vegdist(decostand(dune, method='hellinger'), method='euclid')
disper.dune <- betadisper(dune.dist, dune.env$Management)
anova(disper.dune)
## Analysis of Variance Table
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq F value Pr(>F)
## Groups     3 0.17812 0.059373  1.5729 0.2349
## Residuals 16 0.60397 0.037748
permanova.dune <- adonis2(dune.dist ~ Management, data = dune.env)
permanova.dune
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = dune.dist ~ Management, data = dune.env)
##            Df SumOfSqs      R2      F Pr(>F)   
## Management  3   3.1672 0.29986 2.2842  0.004 **
## Residual   16   7.3950 0.70014                 
## Total      19  10.5621 1.00000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Le test montre que le type d’aménagement a un effet statistiquement significatif sur la composition en espèces végétales des sites.

ggord + ggplot

On va créer un graphique pour présenter la PCA avec les différents types d’aménagement.

Extraire les scores des sites et des espèces, à partir de l’objet créer par prcomp. Cela va nous permettre de mettre éventuellement seulement les espèces désirées sur le graphique.

site.scores<-pca.dune$x 
(site.scores.df <- data.frame(site.scores))
##            PC1         PC2           PC3          PC4          PC5
## 1  -0.42202236 -0.34450749 -3.749271e-01  0.008024023 -0.372799776
## 2  -0.35515437 -0.29665390 -1.142664e-01 -0.042440109 -0.174080679
## 3  -0.08063254 -0.44491000 -3.301696e-02  0.138292966 -0.106457407
## 4  -0.04086906 -0.39718712  4.310306e-02  0.170494078 -0.138027375
## 5  -0.45648854  0.04206103 -8.218385e-02 -0.274349726  0.172225272
## 6  -0.38805443  0.16386896 -5.863489e-02 -0.154889885  0.371372259
## 7  -0.45031849  0.05530804 -3.922469e-02 -0.212042295  0.286577531
## 8   0.29413758 -0.23669691 -3.651072e-02  0.100532017  0.005336318
## 9   0.03585360 -0.39695990  1.557600e-01  0.040285686  0.042171903
## 10 -0.45159580  0.10989979 -1.372343e-01 -0.061864199  0.070651366
## 11 -0.27231082  0.23877731  2.256124e-05  0.444472361  0.062040814
## 12  0.24580687 -0.26730006  5.282509e-01 -0.017884642  0.241545549
## 13  0.22615753 -0.38991817  3.868134e-01 -0.145331015 -0.034885901
## 14  0.56893064  0.20376783 -1.816319e-01 -0.166655594 -0.041816081
## 15  0.65225130  0.23242747 -2.142232e-01  0.005247591  0.096401575
## 16  0.71764913 -0.08613312 -1.564761e-01 -0.182487145  0.003907619
## 17 -0.31641541  0.60652639  1.948300e-01 -0.325939071 -0.393125861
## 18 -0.20046585  0.32084547 -1.151960e-01  0.455690043  0.178975130
## 19  0.00624186  0.58987086  4.687034e-01  0.170712594 -0.246398835
## 20  0.68729916  0.29691351 -2.339574e-01  0.050132321 -0.023613421
##              PC6         PC7          PC8          PC9         PC10
## 1  -0.1947189103  0.14297479  0.194301690  0.112376805  0.032457596
## 2   0.2647456146 -0.17427545 -0.044272021 -0.016720326 -0.143103167
## 3   0.0070264595 -0.07648163 -0.075063932  0.010914507  0.127904407
## 4   0.0420189647 -0.17446434 -0.013889808 -0.162887863  0.140023693
## 5  -0.0492961680 -0.23585887  0.081909418 -0.045412441  0.024452056
## 6  -0.1032250470  0.07513035  0.096314316 -0.011148451  0.095001576
## 7  -0.0582862075  0.08584994  0.002778682  0.046034551 -0.050845625
## 8  -0.0777390258  0.19094461 -0.154446492 -0.153598933 -0.114991844
## 9  -0.1761440885  0.01711273  0.228755610 -0.015221320 -0.171424775
## 10  0.1944992350 -0.05049491 -0.210307980 -0.082508198 -0.058840139
## 11  0.0407933227  0.29970308  0.011721456 -0.056702121  0.077192077
## 12 -0.0074164915 -0.02536257  0.038872424  0.029358985  0.088641851
## 13  0.1359729727  0.14248707 -0.125859881  0.202329767 -0.093788726
## 14  0.5026102903  0.10186304  0.209182998  0.066090412  0.102879773
## 15  0.0009474427  0.02473974  0.061662868 -0.227259133 -0.140607818
## 16 -0.2119652616 -0.02270390 -0.177098525 -0.022669604  0.177063594
## 17 -0.1102816574  0.12212982 -0.151255687 -0.008799111  0.014706687
## 18  0.0362632716 -0.11150681 -0.122134229  0.204232852 -0.018315026
## 19 -0.0359944048 -0.15838135  0.156730498 -0.067611130 -0.008421206
## 20 -0.1998103115 -0.17340534 -0.007901405  0.199200753 -0.079984983
##            PC11         PC12         PC13         PC14         PC15
## 1  -0.054588101  0.034462703 -0.035701587  0.117878451 -0.002213425
## 2   0.019457009  0.121961523 -0.072668906 -0.027268598 -0.114056583
## 3   0.189426459  0.098304109  0.102593714  0.003582509  0.063190553
## 4  -0.139198735 -0.233732456  0.058763244 -0.004413976  0.004692423
## 5   0.048675368 -0.054623058  0.056167712 -0.061130295 -0.067831090
## 6  -0.016507628  0.106475665  0.122588589  0.013653305  0.009980931
## 7  -0.118264835 -0.044205396 -0.022019777  0.012121836  0.002703392
## 8  -0.062119291  0.118226841  0.103731357 -0.096723711 -0.011208506
## 9   0.113994955 -0.061220658 -0.078065568 -0.121604909  0.085266972
## 10 -0.096890649  0.019846981 -0.109878765  0.042312047  0.134363570
## 11 -0.028975224 -0.049496469 -0.072451734 -0.063347298 -0.070411540
## 12  0.042099025  0.020074642 -0.142911945  0.072984250 -0.032959591
## 13 -0.061978791 -0.062186344  0.112696440  0.064723117 -0.016580691
## 14  0.015038111 -0.002557763  0.007135486 -0.066269669  0.032849963
## 15  0.119316583 -0.069303367  0.026613528  0.170846882 -0.023002423
## 16 -0.009598986  0.085246184 -0.110541292 -0.022926853 -0.014930180
## 17  0.148298447 -0.106347535 -0.005376556 -0.035117493 -0.003259164
## 18  0.130118499 -0.034162850  0.015833905  0.027468954  0.003240154
## 19 -0.116765892  0.158724721  0.036465495  0.028445289  0.022972462
## 20 -0.121536323 -0.045487474  0.007026659 -0.055213837 -0.002807226
##             PC16         PC17         PC18         PC19          PC20
## 1   0.0127472182 -0.029903023 -0.024856595 -0.012336708  4.336809e-17
## 2  -0.0446611349  0.005305086  0.022078375  0.031007751  8.673617e-18
## 3  -0.0271552448  0.058753839  0.046914440 -0.030425313  7.979728e-17
## 4  -0.0437358342 -0.010564791 -0.011527629  0.023511185 -4.683753e-17
## 5   0.1097759923 -0.013436149 -0.016784278 -0.042864301 -2.428613e-17
## 6  -0.0306847435 -0.055202850  0.037060951  0.049957778 -1.040834e-17
## 7  -0.0816920274  0.113904666 -0.027202323 -0.012878277 -3.469447e-17
## 8  -0.0336516196 -0.042583349 -0.052191470 -0.028137562  0.000000e+00
## 9   0.0244642383  0.006081289 -0.000449554  0.027768859  7.632783e-17
## 10  0.0459636131 -0.036301400  0.020925899 -0.012264575 -5.724587e-17
## 11  0.0544619654  0.019844250  0.051357700 -0.009062018 -7.806256e-17
## 12 -0.0627871904 -0.055322352 -0.003578525 -0.032477460  5.377643e-17
## 13  0.0862150522  0.009678460  0.014880760  0.014241811  5.724587e-17
## 14 -0.0147797128 -0.003349925 -0.024719508 -0.005808041  2.775558e-17
## 15  0.0061529347  0.016725353  0.012680768  0.001031098  7.285839e-17
## 16  0.0634859914  0.037481480 -0.023923219  0.034616275  4.163336e-17
## 17 -0.0497449145 -0.020609099  0.001575928  0.004009650  8.673617e-18
## 18 -0.0006542089 -0.012827698 -0.067164565  0.018814567 -1.387779e-17
## 19  0.0295200565  0.034403175 -0.015520396  0.002854336  7.459311e-17
## 20 -0.0432404313 -0.022076960  0.060443242 -0.021559055  8.326673e-17
sp.scores<-pca.dune$rotation
(sp.scores.df <- data.frame(sp.scores))
##                    PC1         PC2          PC3          PC4          PC5
## Achimill -0.2248791555  0.07079821 -0.119661537 -0.341760328 -0.101948354
## Agrostol  0.4156189533 -0.22347029  0.034985202 -0.040936713  0.005936338
## Airaprae -0.0351857956  0.19191755  0.202676001 -0.078888514 -0.282855712
## Alopgeni  0.1391150350 -0.39519834  0.324826566  0.002621325  0.006164517
## Anthodor -0.1969422322  0.29608445  0.166870274 -0.353895610 -0.077737618
## Bellpere -0.1152141976 -0.06748736 -0.100435394  0.128114504  0.002135723
## Bromhord -0.1415629580 -0.06594128 -0.086400941 -0.113139773  0.042917985
## Chenalbu  0.0121942554 -0.03216665  0.062905676 -0.030280595 -0.007800813
## Cirsarve -0.0026687294 -0.03968196  0.008489110  0.043021061 -0.037378414
## Comapalu  0.1104463921  0.06035712 -0.107998204 -0.055730602  0.021009915
## Eleopalu  0.3808672608  0.08575416 -0.323855812 -0.126443878  0.029551309
## Elymrepe -0.1461693625 -0.30816535 -0.164224860  0.020107591 -0.302907981
## Empenigr  0.0004910767  0.07100360  0.111218624  0.051899454 -0.080393271
## Hyporadi -0.0560973565  0.24551212  0.242318916  0.072606975 -0.291583945
## Juncarti  0.2486953063 -0.01550957 -0.160744874  0.010897632  0.054225713
## Juncbufo  0.0190979141 -0.15073232  0.312500987 -0.101554267  0.190412207
## Lolipere -0.3377427168 -0.21557574 -0.360101987  0.304110194  0.018829208
## Planlanc -0.2555804426  0.24517118 -0.060083932 -0.078852590  0.300164812
## Poaprat  -0.2818291593 -0.17723703 -0.143688346  0.167481962 -0.102623821
## Poatriv  -0.1321800496 -0.45095346  0.103574859 -0.271062401  0.118784798
## Ranuflam  0.2792578969  0.03027933 -0.141161681 -0.103865168  0.001850673
## Rumeacet -0.1082864119 -0.02990082  0.094139189 -0.241619530  0.430903141
## Sagiproc  0.0410750255 -0.10193982  0.429563593  0.251953139 -0.030424463
## Salirepe  0.0654008835  0.19415423  0.012563845  0.269469132 -0.034009563
## Scorautu -0.0564104250  0.19693367  0.165511278  0.187848971  0.054307016
## Trifprat -0.1004765803  0.03522370 -0.042431386 -0.187403726  0.283990277
## Trifrepe -0.0472757003 -0.00868603  0.092022393  0.090353422  0.285088648
## Vicilath -0.0543676300  0.06549366 -0.040256932  0.226673409  0.078007690
## Bracruta  0.0855581430  0.13397561 -0.006591469  0.358036940  0.449131491
## Callcusp  0.2051900747  0.07088792 -0.181340715 -0.128624446 -0.029851304
##                   PC6           PC7          PC8         PC9         PC10
## Achimill -0.001433447 -4.245252e-05 -0.139275612 -0.03103964 -0.144243870
## Agrostol  0.046045557 -3.525837e-02 -0.078558371 -0.10461922  0.210297407
## Airaprae -0.089716769 -1.169179e-02 -0.018998063 -0.09251269  0.013279470
## Alopgeni -0.051261763  1.777523e-03 -0.380074914 -0.02631089  0.158579434
## Anthodor -0.112317801 -1.080043e-01 -0.106432326 -0.23023176  0.031200744
## Bellpere  0.213317439 -4.816241e-01 -0.271368034 -0.03193991  0.036080560
## Bromhord  0.223494039 -3.649116e-01 -0.185180260 -0.27427582 -0.154746474
## Chenalbu  0.041261458  6.203707e-02 -0.064291687  0.13439007 -0.078827421
## Cirsarve  0.015441982 -9.199161e-02 -0.008592696 -0.13102732  0.142525942
## Comapalu  0.253410826  9.179244e-02  0.230556545 -0.18290679 -0.056799466
## Eleopalu  0.008558863  1.001312e-01 -0.072559072 -0.25624816 -0.007112988
## Elymrepe -0.134141352 -3.202739e-01  0.473903145 -0.06448830 -0.019325678
## Empenigr -0.015937423 -1.006169e-01  0.116818638 -0.06552646 -0.010327426
## Hyporadi -0.077618969  1.398465e-01  0.031234356 -0.16995428  0.102773191
## Juncarti -0.373541740  1.366486e-02 -0.035830116 -0.26948199 -0.557137446
## Juncbufo -0.050382648  1.472281e-01  0.136186952  0.29199405 -0.302167976
## Lolipere -0.077527961  2.633129e-01  0.057940579 -0.13474471  0.104356918
## Planlanc -0.061099106  1.501356e-01 -0.257949680  0.08738485  0.137031409
## Poaprat  -0.066342051  2.883555e-01 -0.173483287  0.08478524 -0.165446302
## Poatriv  -0.032482359 -9.052647e-02 -0.169807161  0.02929114 -0.104717464
## Ranuflam  0.065379870  1.165436e-01 -0.097973268  0.12783117 -0.220685143
## Rumeacet -0.190843551 -8.174166e-02  0.357868885 -0.01191780  0.056890788
## Sagiproc -0.020650816  1.011335e-01  0.130867029 -0.25731110  0.034324439
## Salirepe -0.138332766 -3.903741e-01  0.014296262  0.48475738 -0.197219009
## Scorautu  0.245056856 -6.067166e-02 -0.046290901 -0.05730164 -0.333608160
## Trifprat -0.099328659 -1.856270e-02  0.144877989 -0.01182233  0.118607508
## Trifrepe  0.603878190  5.671622e-02  0.209187168 -0.13289241 -0.132042485
## Vicilath  0.081648423  1.144655e-01 -0.154486392  0.04787343  0.032832865
## Bracruta -0.332896166 -2.195302e-01 -0.081168769 -0.21669754  0.145243528
## Callcusp  0.137925615 -4.803121e-02  0.086693108  0.31331767  0.340411008
##                   PC11         PC12         PC13          PC14         PC15
## Achimill -1.255743e-02  0.046215103 -0.256723257  0.1313444344 -0.057272752
## Agrostol  2.532151e-02 -0.420175415  0.252885419 -0.0253551477  0.347506784
## Airaprae  9.391228e-02  0.060848890  0.077436256 -0.0404887435  0.106565696
## Alopgeni  2.913598e-01  0.340305846 -0.063431110 -0.2035060907 -0.108231449
## Anthodor -5.605106e-02  0.047468359  0.162324336 -0.0767398112  0.536035731
## Bellpere  2.277214e-01 -0.066340892  0.071266730 -0.0407483438 -0.002383024
## Bromhord -3.973954e-01 -0.221153049 -0.277170883 -0.0725695655 -0.125790738
## Chenalbu -5.683765e-02 -0.062470024  0.161943974  0.1147861266 -0.051639846
## Cirsarve -1.545945e-01 -0.284355335  0.102264776 -0.0094803775  0.017698853
## Comapalu  2.082232e-01 -0.122194920  0.081787344  0.3183693811  0.048304841
## Eleopalu -3.295834e-02  0.171174153 -0.031189575 -0.0927636652 -0.164919807
## Elymrepe  2.981940e-01 -0.141465524 -0.013771024 -0.1796833253 -0.059413586
## Empenigr -1.562428e-01  0.232654972  0.076458863  0.0736090906  0.104395306
## Hyporadi  6.753349e-05  0.072357556 -0.044834636 -0.1755997380 -0.171165046
## Juncarti  6.361967e-02  0.016314991 -0.103042130 -0.3378149280  0.160162220
## Juncbufo  2.254486e-02 -0.235107742 -0.357842523  0.0954648841  0.192810046
## Lolipere -2.665899e-01  0.181612509  0.033499240 -0.0317211176  0.151252550
## Planlanc  1.713019e-01 -0.346381680  0.025105687 -0.2449153380 -0.091817541
## Poaprat   8.882030e-02 -0.157511234  0.260442751 -0.2080073931  0.200952629
## Poatriv  -2.214968e-01  0.194808307  0.181781411  0.0678421701  0.010998379
## Ranuflam -1.877691e-01 -0.031168655  0.298475438  0.0007910441 -0.153366894
## Rumeacet  7.013359e-02 -0.009512001  0.043474284 -0.2219769387 -0.145511466
## Sagiproc -4.214871e-01 -0.209789606 -0.013130855 -0.1631780935 -0.143417999
## Salirepe -2.200017e-01  0.113806094  0.160506778 -0.0424753661  0.127010328
## Scorautu  1.534068e-01 -0.042566093  0.244233847 -0.2416944201 -0.275228477
## Trifprat -1.120780e-01  0.073288375  0.385957343 -0.0618061194 -0.193278565
## Trifrepe  7.257752e-02  0.289166907  0.008895578 -0.1890818479  0.325637771
## Vicilath  1.591913e-02 -0.091882888 -0.262687703 -0.0417491662  0.062814819
## Bracruta  6.889939e-02  0.075253778 -0.096255783  0.2311884646  0.117003891
## Callcusp -1.820802e-01  0.060639205 -0.233039917 -0.5210466770  0.143774250
##                  PC16        PC17        PC18        PC19        PC20
## Achimill -0.267107657 -0.41210202  0.22489767  0.09626702 -0.07897943
## Agrostol -0.199229193  0.04314512  0.29537634 -0.10446748  0.06311722
## Airaprae -0.172178271  0.10130470 -0.19670085  0.19247304  0.18122527
## Alopgeni -0.192503646  0.01423069 -0.03731558 -0.11972890  0.14908814
## Anthodor  0.124188453 -0.05692455 -0.01459231 -0.33620497 -0.11244598
## Bellpere  0.118498260 -0.06025397 -0.15845467 -0.02686258 -0.42878365
## Bromhord -0.108289301  0.33193794  0.02391451 -0.01840574  0.14495526
## Chenalbu  0.287726031  0.05372402  0.11981393  0.20287560 -0.02388393
## Cirsarve -0.176766183 -0.07102134 -0.11240559  0.40560643  0.32331465
## Comapalu -0.047010757  0.12643349 -0.15710049 -0.11232080 -0.19504909
## Eleopalu  0.036796973  0.11125275 -0.49725476 -0.12181591 -0.03419597
## Elymrepe  0.255508256  0.03734245  0.05576828 -0.11750010  0.15496839
## Empenigr  0.143748846  0.27864541 -0.18233754  0.05932814  0.20581585
## Hyporadi  0.140080500  0.35880767  0.33217606  0.02822715 -0.13068216
## Juncarti  0.052524852 -0.06945109  0.11264250  0.22398356 -0.02128059
## Juncbufo -0.114034976  0.36869437 -0.13618729 -0.08144971 -0.16453618
## Lolipere -0.272395734  0.13416285  0.10854839 -0.33124013  0.02818100
## Planlanc  0.174097492  0.08161027 -0.18223409 -0.10838544  0.31479145
## Poaprat  -0.036402138  0.05182841 -0.23735935  0.22250549 -0.20920638
## Poatriv   0.295850357  0.12398510  0.01209593  0.03105173 -0.03980571
## Ranuflam  0.217490753 -0.05986007  0.20438397 -0.27663575  0.14252283
## Rumeacet -0.104640190 -0.15313263 -0.04734271 -0.17877450  0.07037098
## Sagiproc  0.115654700 -0.39453594 -0.23802349 -0.10195268 -0.24920446
## Salirepe -0.161050639 -0.07780178 -0.13606538 -0.12265427  0.13147778
## Scorautu -0.280606173  0.13110453  0.18988616 -0.18533745 -0.08861247
## Trifprat -0.086183943  0.15164056  0.10448329  0.32730883 -0.23536861
## Trifrepe  0.084685940 -0.12711462  0.04985696  0.19051210  0.20983191
## Vicilath  0.392993308 -0.09705062  0.14360459 -0.04214081  0.11269105
## Bracruta  0.058168027  0.09580990  0.15422886  0.09217431 -0.10719603
## Callcusp -0.006585697  0.09775620  0.06929603  0.11123932 -0.28379314

Ouverture des librairies nécessaires à la création du graphique. Si problème avec installation de la librairie ggord : https://rdrr.io/github/fawda123/ggord/

library(ggord)
library(ggplot2)
library(ggrepel) # pour éviter overlap du texte

On va d’abord commencer le graphique avec la fonctions ggord du package du même nom.

Le graphique créé par la fonction ggord est un objet ggplot. Cet objet (ici appelé pca.plot) peut être personnalisé à même la fonction ggord, et peut aussi servir de base sur laquelle ajouter des couches avec ggplot.

pca.plot <- ggord(pca.dune, 
                  xlims=c(-1,1), ylims=c(-1,1)) # xlims et ylims définissent les limites dans l'affichage des axe (les valeurs entre parenthèses correspondent à l'échelle numérique des axes). On ajuste au fur et à mesure, au besoin.
pca.plot # À la fin de chaque chunk, l'objet graphique est appelé pour qu'on visualise le rendu à chaque fois.

On remarque que la fonctionggord inclut par défaut la proportion de la variance expliquée par les axes. C’est super!

Par défaut, la fonction ggord va mettre en graphique les deux premiers axes de l’ordination. Mais si on s’intéresse à différents axes, on peut spécifier lesquels mettre en graphique.

pca.plot <- ggord(pca.dune, 
                  xlims=c(-1,1), ylims=c(-1,1),
                  axes=c(1,3)) # Ici, on met en graphique les axes 1 et 3.
pca.plot

On revient aux deux premiers axes. Faisons quelques petits changements. Les changements sont marqués d’un “#” à côté de la ligne.

pca.plot <- ggord(pca.dune, 
                  xlims=c(-1.2,1.4), ylims=c(-0.75,0.9), # ajustement des limites
                  axes=c(1,2), # retour aux axes 1 et 2
                  size=3, # diminution de la taille des points (sites)
                  arrow=NULL,# enlever les flèches (qui viennent par défaut) pour les espèces
                  labcol="forestgreen") # changer la couleur du texte pour les espèces
pca.plot

Faisons d’autres changements!

pca.plot <- ggord(pca.dune, 
                  xlims=c(-1.2,1.4), ylims=c(-0.75,0.9), 
                  axes=c(1,2), 
                  size=3, 
                  obslab=TRUE, # pour passer de points pour les sites (par défaut) à texte (nom des sites)
                  arrow=NULL,
                  labcol="forestgreen")
pca.plot

Il y a beaucoup d’espèces affichées au milieu et c’est illisible. En plus, les espèces au milieu ne contribuent pas beaucoup à la distinction compositionnelle des sites. On va donc garder seulement les quelques espèces dont les scores ont la valeur absolue la plus élevée, pour chacun des deux axes, qu’on met en graphique (ici PC1 et PC2).

D’abord, on va créer un nouveau dataframe contenant seulement ces espèces, en utilisant quelques fonctions de la librairie dplyr.

library(dplyr)
A <- top_n(sp.scores.df, 3, PC1) # sélectionne les 3 espèces dont le score est le plus élevé sur l'axe 1.
B <- top_n(sp.scores.df, -3, PC1) # sélectionne les 3 espèces dont le score est le plus bas sur l'axe 1.
C <- top_n(sp.scores.df, 3, PC2) # sélectionne les 3 espèces dont le score est le plus élevé sur l'axe 2.
D <- top_n(sp.scores.df, -3, PC2) # sélectionne les 3 espèces dont le score est le plus bas sur l'axe 2.
sp.scores.skim <- bind_rows(A,B,C,D) # merge les sub dataframe A, B, C et D.
(sp.scores.skim <- distinct(sp.scores.skim)) # on ne garde que les rangées (espèces) qui sont uniques (pour éviter qu'une même espèce s'y retrouve en copies)
##                  PC1         PC2         PC3          PC4          PC5
## Agrostol  0.41561895 -0.22347029  0.03498520 -0.040936713  0.005936338
## Eleopalu  0.38086726  0.08575416 -0.32385581 -0.126443878  0.029551309
## Ranuflam  0.27925790  0.03027933 -0.14116168 -0.103865168  0.001850673
## Lolipere -0.33774272 -0.21557574 -0.36010199  0.304110194  0.018829208
## Planlanc -0.25558044  0.24517118 -0.06008393 -0.078852590  0.300164812
## Poaprat  -0.28182916 -0.17723703 -0.14368835  0.167481962 -0.102623821
## Anthodor -0.19694223  0.29608445  0.16687027 -0.353895610 -0.077737618
## Hyporadi -0.05609736  0.24551212  0.24231892  0.072606975 -0.291583945
## Alopgeni  0.13911504 -0.39519834  0.32482657  0.002621325  0.006164517
## Elymrepe -0.14616936 -0.30816535 -0.16422486  0.020107591 -0.302907981
## Poatriv  -0.13218005 -0.45095346  0.10357486 -0.271062401  0.118784798
##                   PC6          PC7         PC8         PC9         PC10
## Agrostol  0.046045557 -0.035258372 -0.07855837 -0.10461922  0.210297407
## Eleopalu  0.008558863  0.100131193 -0.07255907 -0.25624816 -0.007112988
## Ranuflam  0.065379870  0.116543559 -0.09797327  0.12783117 -0.220685143
## Lolipere -0.077527961  0.263312916  0.05794058 -0.13474471  0.104356918
## Planlanc -0.061099106  0.150135587 -0.25794968  0.08738485  0.137031409
## Poaprat  -0.066342051  0.288355460 -0.17348329  0.08478524 -0.165446302
## Anthodor -0.112317801 -0.108004332 -0.10643233 -0.23023176  0.031200744
## Hyporadi -0.077618969  0.139846549  0.03123436 -0.16995428  0.102773191
## Alopgeni -0.051261763  0.001777523 -0.38007491 -0.02631089  0.158579434
## Elymrepe -0.134141352 -0.320273944  0.47390315 -0.06448830 -0.019325678
## Poatriv  -0.032482359 -0.090526465 -0.16980716  0.02929114 -0.104717464
##                   PC11        PC12        PC13          PC14        PC15
## Agrostol  2.532151e-02 -0.42017542  0.25288542 -0.0253551477  0.34750678
## Eleopalu -3.295834e-02  0.17117415 -0.03118957 -0.0927636652 -0.16491981
## Ranuflam -1.877691e-01 -0.03116866  0.29847544  0.0007910441 -0.15336689
## Lolipere -2.665899e-01  0.18161251  0.03349924 -0.0317211176  0.15125255
## Planlanc  1.713019e-01 -0.34638168  0.02510569 -0.2449153380 -0.09181754
## Poaprat   8.882030e-02 -0.15751123  0.26044275 -0.2080073931  0.20095263
## Anthodor -5.605106e-02  0.04746836  0.16232434 -0.0767398112  0.53603573
## Hyporadi  6.753349e-05  0.07235756 -0.04483464 -0.1755997380 -0.17116505
## Alopgeni  2.913598e-01  0.34030585 -0.06343111 -0.2035060907 -0.10823145
## Elymrepe  2.981940e-01 -0.14146552 -0.01377102 -0.1796833253 -0.05941359
## Poatriv  -2.214968e-01  0.19480831  0.18178141  0.0678421701  0.01099838
##                 PC16        PC17        PC18        PC19        PC20
## Agrostol -0.19922919  0.04314512  0.29537634 -0.10446748  0.06311722
## Eleopalu  0.03679697  0.11125275 -0.49725476 -0.12181591 -0.03419597
## Ranuflam  0.21749075 -0.05986007  0.20438397 -0.27663575  0.14252283
## Lolipere -0.27239573  0.13416285  0.10854839 -0.33124013  0.02818100
## Planlanc  0.17409749  0.08161027 -0.18223409 -0.10838544  0.31479145
## Poaprat  -0.03640214  0.05182841 -0.23735935  0.22250549 -0.20920638
## Anthodor  0.12418845 -0.05692455 -0.01459231 -0.33620497 -0.11244598
## Hyporadi  0.14008050  0.35880767  0.33217606  0.02822715 -0.13068216
## Alopgeni -0.19250365  0.01423069 -0.03731558 -0.11972890  0.14908814
## Elymrepe  0.25550826  0.03734245  0.05576828 -0.11750010  0.15496839
## Poatriv   0.29585036  0.12398510  0.01209593  0.03105173 -0.03980571

Mise en graphique de ces espèces seulement

pca.plot <- ggord(pca.dune, 
                  xlims=c(-1.2,1.4), ylims=c(-0.75,0.9), 
                  axes=c(1,2), 
                  size=3, 
                  obslab=TRUE,
                  arrow=NULL,
                  txt=NULL, # enlever le texte (par défaut) pour les espèces
                  addpts=sp.scores.skim, # Ajouter les espèces aux plus hauts et faibles scores
                  addsize=3, # taille des points des espèces
                  addcol="forestgreen") # couleur des espèces
                  # on enlève l'argument labcol
pca.plot

pca.plot <- ggord(pca.dune, 
                  xlims=c(-1.2,1.4), ylims=c(-0.75,0.9), 
                  axes=c(1,2), 
                  size=3, 
                  obslab=TRUE,
                  arrow=NULL,
                  txt=NULL, 
                  addpts=sp.scores.skim, 
                  addsize=4, # augmenter la taille du texte des espèces
                  ptslab=TRUE, # changer de points à texte pour les espèces
                  addcol="forestgreen")
pca.plot

On va ajouter les résultats d’un envfit. D’abord, on fait le test, en s’assurant qu’on calcule la régression sur les deux axes qu’on a choisit de mettre en graphique (par défaut, 1 et 2).

(pca.envfit <- envfit(pca.dune, dune.env))
## 
## ***VECTORS
## 
##        PC1     PC2     r2 Pr(>r)  
## A1 0.97222 0.23409 0.3622  0.019 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
## 
## ***FACTORS:
## 
## Centroids:
##                  PC1     PC2
## Moisture1    -0.3635  0.0257
## Moisture2    -0.2224 -0.0314
## Moisture4     0.1408 -0.3321
## Moisture5     0.4504  0.0872
## ManagementBF -0.3597  0.0173
## ManagementHF -0.1930 -0.0745
## ManagementNM  0.2330  0.3751
## ManagementSF  0.1077 -0.3217
## UseHayfield  -0.0994  0.2242
## UseHaypastu  -0.0203 -0.2180
## UsePasture    0.1716  0.0350
## Manure0       0.2330  0.3751
## Manure1      -0.2294 -0.0161
## Manure2      -0.2385 -0.0895
## Manure3       0.1969 -0.1644
## Manure4      -0.1812 -0.3955
## 
## Goodness of fit:
##                r2 Pr(>r)    
## Moisture   0.5366  0.001 ***
## Management 0.4614  0.003 ** 
## Use        0.1794  0.134    
## Manure     0.4531  0.011 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999

Pour incorporer les vecteurs résultant de ce test, il faut d’abord les extraire et les mettre dans un dataframe.

vect <- data.frame(pca.envfit[["vectors"]][["arrows"]])
vect$variable.env<-rownames(vect) # il mettre le rowname(la variable environnementale en colonne)
vect
##          PC1       PC2 variable.env
## A1 0.9722159 0.2340862           A1
pca.plot <- ggord(pca.dune, 
                  xlims=c(-1.2,1.4), ylims=c(-0.75,0.9), 
                  axes=c(1,2), 
                  size=3, 
                  obslab=TRUE,
                  arrow=NULL,
                  txt=NULL, 
                  addpts=sp.scores.skim, 
                  addsize=4, 
                  ptslab=TRUE, 
                  addcol="forestgreen")
pca.plot + 
  coord_fixed() + # ajout des vecteurs de envfit
  geom_segment(data = vect,
               aes(x = 0, xend = PC1, y = 0, yend = PC2),
               arrow = arrow(length = unit(0.5, "cm")),
               colour = "darkgrey") +
  geom_text_repel(data = vect, 
                  aes(x = PC1, y = PC2, label = variable.env, size = 3), 
                  show.legend = FALSE)
## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.

Maintenant, on va ajouter des ellipses autour des sites appartenant à des groupes similaires. La fonction ggord est à même de calculer ces ellipses avec intervalle de confiance spécifiée. On va aussi revenir à des symboles au lieu du nom des sites et changer la couleur du texte des espèces.

pca.plot <- ggord(pca.dune, 
                  xlims=c(-1.2,1.4), ylims=c(-0.75,0.9),
                  axes=c(1,2),
                  grp_in=dune.env$Management, # ajout des ellipses en fonction du type de Management
                  grp_title = "Management",
                  ellipse_pro=0.95, # valeur de confiance pour les ellipses
                  alpha_el=0.3, # transparence des ellipses
                  size=3, 
                  obslab=FALSE, # changement de texte à points pour sites
                  arrow=NULL,
                  txt=NULL, 
                  addpts=sp.scores.skim, 
                  addsize=4, 
                  ptslab=TRUE,
                  addcol="black", # couleur des espèces
                  alpha = 1)
pca.plot + 
    scale_shape_manual('Groups', values = c(15,16,17,18)) + # Symboles pour les sites en fonction de leur groupe. Assurer vous d'avoir un nombre de symboles correspondant au nombre de groupes.
    coord_fixed() +
    geom_segment(data = vect,
               aes(x = 0, xend = PC1, y = 0, yend = PC2),
               arrow = arrow(length = unit(0.5, "cm")), 
               colour = "darkgrey") +
    geom_text_repel(data = vect, 
               aes(x = PC1, y = PC2, label = variable.env, size = 3), 
               show.legend = FALSE)
## Scale for shape is already present.
## Adding another scale for shape, which will replace the existing scale.
## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.

Voici une légende indiquant quels symboles sont associées aux différentes valeurs numériques de la couche scale_shape_manual.

Changer les couleurs des groupes (types d’aménagement). On va extraire les codes des couleurs d’une palette de RColorBrewer avec la fonction brewer.pal

library(RColorBrewer)
pca.colors <- brewer.pal(n = 4, name = "Dark2") # Codes de couleur héxadécimaux

On va incorporer ces codes de couleurs dans un argument. On déplace la légende vers le haut.

pca.plot <- ggord(pca.dune, 
                  xlims=c(-1.2,1.4), ylims=c(-0.75,0.9),
                  axes=c(1,2),
                  grp_in=dune.env$Management, 
                  grp_title = "Management",
                  cols=pca.colors, # ajout des couleurs choisies
                  ellipse_pro=0.95,
                  alpha_el=0.3,
                  size=3, 
                  obslab=FALSE, 
                  arrow=NULL,
                  txt=NULL, 
                  addpts=sp.scores.skim, 
                  addsize=4, 
                  ptslab=TRUE,
                  addcol="black",
                  alpha = 1)
pca.plot <- pca.plot + # ici j'enregistre l'ensemble des arguments précédents et suivants dans mon objet pca.plot
    scale_shape_manual('Groups', values = c(15,16,17,18)) + 
    coord_fixed() +
    geom_segment(data = vect,
               aes(x = 0, xend = PC1, y = 0, yend = PC2),
               arrow = arrow(length = unit(0.5, "cm")), 
               colour = "darkgrey") +
    geom_text_repel(data = vect, 
                    aes(x = PC1, y = PC2, label = variable.env, size = 3), 
                    show.legend = FALSE) +
    theme(legend.position = 'top') # mettre la légende en haut
## Scale for shape is already present.
## Adding another scale for shape, which will replace the existing scale.
## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.
pca.plot

Créer un thème custom pour ggplot

Allez explorer les options (liste des arguments dans le fichier help, lorsque vous rentrez ?theme dans la console de R).

super_theme <- theme(
        panel.background = element_blank(),
        panel.grid = element_blank(), # grillage, ici aucun, mais on aurait pu mettre par ex.  element_line("black"),
        axis.line = element_line("black"),
        text = element_text(size = 12),
        axis.text = element_text(size = 10, colour = "gray25"), # taille des valeurs des axes
        axis.title = element_text(size = 14, colour = "gray25"), # taille des titres des axes
        legend.title = element_text(size = 14),
        legend.text = element_text(size = 14),
        legend.key = element_blank())

Ajouter ce thème au graphique.

pca.plot + super_theme

On va rajouter du texte!

pca.plot + 
    annotate("text", x = 1.1, y = 0.7, label = "PERMANOVA\nR2=30%\n p=0.007") +
    super_theme

Sauvegarder

Satisfait? On sauvegarde l’image en png. Par défaut, ggplot enregistre le dernier graphique créé. Ici on copie-colle le code du graphique ci haut, puis sauvegarde ce graphique avec la fonction ggsave dans un dossier qu’on a créé (figures/).

pca.plot + 
    annotate("text", x = 1.1, y = 0.7, label = "PERMANOVA\nR2=30%\n p=0.007") +
    super_theme

ggsave("pca.png", path="./figures/")
## Saving 7 x 5 in image

Cet atelier n’est pas fait pour vous montrer à utiliser une esthétique particulière, mais bien pour vous familiariser avec les options (presque infinies…) permettant de personnaliser un graphique d’ordination. Allez voir des figures d’ordination dans la littérature. Essayez de reproduire l’esthétique qui vous semble la plus appropriée… ou, tout simplement, amusez-vous les artistes!!

CA

Préparation des données

Faire une CA avec la fonction cca.

ca.dune <- cca(dune)
ca.summary <- summary(ca.dune)
summary(ca.dune)
## 
## Call:
## cca(X = dune) 
## 
## Partitioning of scaled Chi-square:
##               Inertia Proportion
## Total           2.115          1
## Unconstrained   2.115          1
## 
## Eigenvalues, and their contribution to the scaled Chi-square 
## 
## Importance of components:
##                          CA1    CA2    CA3     CA4     CA5     CA6     CA7
## Eigenvalue            0.5360 0.4001 0.2598 0.17598 0.14476 0.10791 0.09247
## Proportion Explained  0.2534 0.1892 0.1228 0.08319 0.06844 0.05102 0.04372
## Cumulative Proportion 0.2534 0.4426 0.5654 0.64858 0.71702 0.76804 0.81175
##                           CA8     CA9    CA10    CA11    CA12    CA13     CA14
## Eigenvalue            0.08091 0.07332 0.05630 0.04826 0.04125 0.03523 0.020529
## Proportion Explained  0.03825 0.03466 0.02661 0.02282 0.01950 0.01665 0.009705
## Cumulative Proportion 0.85000 0.88467 0.91128 0.93410 0.95360 0.97025 0.979955
##                           CA15     CA16     CA17     CA18     CA19
## Eigenvalue            0.014911 0.009074 0.007938 0.007002 0.003477
## Proportion Explained  0.007049 0.004290 0.003753 0.003310 0.001644
## Cumulative Proportion 0.987004 0.991293 0.995046 0.998356 1.000000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## 
## 
## Species scores
## 
##               CA1      CA2      CA3       CA4       CA5      CA6
## Achimill -0.90859  0.08461 -0.58636 -0.008919 -0.660183 -0.18877
## Agrostol  0.93378 -0.20651  0.28165  0.024293 -0.139326 -0.02256
## Airaprae -1.00434  3.06749  1.33773  0.194305 -1.081813 -0.53699
## Alopgeni  0.40088 -0.61839  0.85013  0.346740  0.016509  0.10169
## Anthodor -0.96676  1.08361 -0.17188  0.459788 -0.607533 -0.30425
## Bellpere -0.50018 -0.35503 -0.15239 -0.704153 -0.058546  0.07308
## Bromhord -0.65762 -0.40634 -0.30685 -0.496751 -0.561358  0.07004
## Chenalbu  0.42445 -0.84402  1.59029  1.248755 -0.207480  0.87566
## Cirsarve -0.05647 -0.76398  0.91793 -1.175919 -0.384024 -0.13985
## Comapalu  1.91690  0.52150 -1.18215 -0.021738 -1.359988  1.31207
## Eleopalu  1.76383  0.34562 -0.57336 -0.002976 -0.332396 -0.14688
## Elymrepe -0.37074 -0.74148  0.26238 -0.566308 -0.270122 -0.72624
## Empenigr -0.69027  3.26420  1.95716 -0.176936 -0.073518 -0.16083
## Hyporadi -0.85408  2.52821  1.13951 -0.175115 -0.311874  0.11177
## Juncarti  1.27580  0.09963 -0.09320  0.005536  0.289410 -0.78247
## Juncbufo  0.08157 -0.68074  1.00545  1.078390  0.268360  0.24168
## Lolipere -0.50272 -0.35955 -0.21821 -0.474727  0.101494 -0.01594
## Planlanc -0.84058  0.24886 -0.78066  0.371149  0.271377  0.11989
## Poaprat  -0.38919 -0.32999 -0.02015 -0.358371  0.079296 -0.05165
## Poatriv  -0.18185 -0.53997  0.23388  0.178834 -0.155342 -0.07584
## Ranuflam  1.55886  0.30700 -0.29765  0.046974 -0.008747 -0.14744
## Rumeacet -0.65289 -0.25525 -0.59728  1.160164  0.255849 -0.32730
## Sagiproc  0.00364  0.01719  1.11570  0.066981  0.186654  0.32463
## Salirepe  0.61035  1.54868  0.04970 -0.607136  1.429729 -0.55183
## Scorautu -0.19566  0.38884  0.03975 -0.130392  0.141232  0.23717
## Trifprat -0.88116 -0.09792 -1.18172  1.282429  0.325706 -0.33388
## Trifrepe -0.07666 -0.02032 -0.20594  0.026462 -0.186748  0.53957
## Vicilath -0.61893  0.37140 -0.46057 -1.000375  1.162652  1.44971
## Bracruta  0.18222  0.26477 -0.16606  0.064009  0.576334  0.07741
## Callcusp  1.95199  0.56743 -0.85948 -0.098969 -0.556737  0.23282
## 
## 
## Site scores (weighted averages of species scores)
## 
##         CA1        CA2      CA3      CA4      CA5      CA6
## 1  -0.81167 -1.0826714 -0.14479 -2.10665 -0.39287 -1.83462
## 2  -0.63268 -0.6958357 -0.09708 -1.18695 -0.97686  0.06575
## 3  -0.10148 -0.9128732  0.68815 -0.68137 -0.08709 -0.28678
## 4  -0.05647 -0.7639784  0.91793 -1.17592 -0.38402 -0.13985
## 5  -0.95293 -0.1846015 -0.95609  0.86853 -0.34552 -0.98333
## 6  -0.85633 -0.0005408 -1.39735  1.59909  0.65494 -0.19386
## 7  -0.87149 -0.2547040 -0.86830  0.90468  0.17385 -0.03446
## 8   0.76268 -0.2968459  0.35648 -0.10772  0.17507 -0.36444
## 9   0.09693 -0.7864314  0.86492  0.40090  0.28704 -1.02783
## 10 -0.87885 -0.0353136 -0.82987 -0.68053 -0.75438  0.81070
## 11 -0.64223  0.4440332 -0.17371 -1.09684  1.37462  2.00626
## 12  0.28557 -0.6656161  1.64423  1.71496  0.65381  1.17376
## 13  0.42445 -0.8440195  1.59029  1.24876 -0.20748  0.87566
## 14  1.91996  0.5351062 -1.39863 -0.08575 -2.21317  2.43044
## 15  1.91384  0.5079036 -0.96567  0.04227 -0.50681  0.19370
## 16  2.00229  0.1090627 -0.33414  0.33760 -0.50097 -0.76159
## 17 -1.47545  2.7724102  0.40859  0.75117 -2.59425 -1.10122
## 18 -0.31241  0.6328355 -0.66501 -1.12728  2.65575  0.97565
## 19 -0.69027  3.2642026  1.95716 -0.17694 -0.07352 -0.16083
## 20  1.94438  1.0688809 -0.66595 -0.55317  1.59606 -1.70292

Extraction de la proportion de la variance expliquée par les deux premiers axes (pour s’y référer plus rapidement).

(prop.expl.ca1 <- ca.summary$cont$importance[2, "CA1"])
## [1] 0.2533987
(prop.expl.ca2 <- ca.summary$cont$importance[2, "CA2"])
## [1] 0.1891696

envfit

(ca.envfit <- envfit(ca.dune, dune.env))
## 
## ***VECTORS
## 
##         CA1      CA2     r2 Pr(>r)  
## A1 0.998160 0.060614 0.3104   0.05 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
## 
## ***FACTORS:
## 
## Centroids:
##                  CA1     CA2
## Moisture1    -0.7484 -0.1423
## Moisture2    -0.4652 -0.2156
## Moisture4     0.1827 -0.7315
## Moisture5     1.1143  0.5708
## ManagementBF -0.7258 -0.1413
## ManagementHF -0.3867 -0.2960
## ManagementNM  0.6517  1.4405
## ManagementSF  0.3376 -0.6761
## UseHayfield  -0.2861  0.6488
## UseHaypastu  -0.0735 -0.5602
## UsePasture    0.5163  0.0508
## Manure0       0.6517  1.4405
## Manure1      -0.4639 -0.1738
## Manure2      -0.5872 -0.3600
## Manure3       0.5187 -0.3172
## Manure4      -0.2059 -0.8775
## 
## Goodness of fit:
##                r2 Pr(>r)    
## Moisture   0.4113  0.006 ** 
## Management 0.4441  0.001 ***
## Use        0.1845  0.092 .  
## Manure     0.4552  0.005 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999

Pour incorporer les vecteurs résultant de ce test, il faut d’abord les extraire et les mettre dans un dataframe.

vect.ca.env <- data.frame(ca.envfit[["vectors"]][["arrows"]])
vect.ca.env$variable.env <- rownames(vect.ca.env)
vect.ca.env
##          CA1        CA2 variable.env
## A1 0.9981613 0.06061359           A1

PERMANOVA. On va tester si la composition végétale diffère entre différents type d’aménagement (facteur; variable quantitative) des sites.

dune.dist.chi2 <- vegdist(dune, method='chisq')
disper.dune.chi2 <- betadisper(dune.dist.chi2, dune.env$Management)
anova(disper.dune.chi2)
## Analysis of Variance Table
## 
## Response: Distances
##           Df Sum Sq Mean Sq F value   Pr(>F)   
## Groups     3 3.0277  1.0092  5.9124 0.006493 **
## Residuals 16 2.7311  0.1707                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
permanova.dune.ca <- adonis2(dune.dist.chi2 ~ Management, data = dune.env)
permanova.dune.ca
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = dune.dist.chi2 ~ Management, data = dune.env)
##            Df SumOfSqs      R2      F Pr(>F)   
## Management  3   12.399 0.25508 1.8263  0.006 **
## Residual   16   36.208 0.74492                 
## Total      19   48.606 1.00000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ici, le résultat du test de dispersion indique que la différence dans la composition végétale entre les types d’aménagement qui est calculée par la PERMANOVA pourrait être du à la différence dans la variance interne en composition des groupes (les groupes ont une dispersion hétérogène).

Extraire les scores de CA avec la fonction scores. Pour spécifier à ggplot quels sont les “points” ou le “texte” à mettre en graphique.

site.scores.dune.ca <- scores(ca.dune, display="sites", choices=c(1,2)) # ici on garde que les scores des axes 1 et 2 mais on aurait pu extraire tous les scores aussi.
site.scores.dune.ca <- data.frame(site.scores.dune.ca)
Management <- dune.env$Management
(site.scores.ca <- cbind(site.scores.dune.ca, Management)) # incorporer variable Management dans les scores (pour aider plus loin dans ggplot à changer l'aspect visuel des points des sites en fonction de leur Management)
##            CA1           CA2 Management
## 1  -0.81167372 -1.0826713631         SF
## 2  -0.63267723 -0.6958357018         BF
## 3  -0.10147537 -0.9128731535         SF
## 4  -0.05647118 -0.7639783770         SF
## 5  -0.95292777 -0.1846015038         HF
## 6  -0.85632779 -0.0005407796         HF
## 7  -0.87149118 -0.2547039586         HF
## 8   0.76268386 -0.2968458686         HF
## 9   0.09693281 -0.7864314398         HF
## 10 -0.87885097 -0.0353135691         BF
## 11 -0.64223343  0.4440332251         BF
## 12  0.28557236 -0.6656161454         SF
## 13  0.42445102 -0.8440194551         SF
## 14  1.91996119  0.5351062139         NM
## 15  1.91384431  0.5079035784         NM
## 16  2.00229065  0.1090626959         SF
## 17 -1.47545057  2.7724102369         NM
## 18 -0.31241136  0.6328354520         NM
## 19 -0.69026877  3.2642026273         NM
## 20  1.94438046  1.0688808879         NM
sp.scores.dune.ca <- scores(ca.dune, display="species", choices=c(1,2)) # idem scores axes 1 et 2 seulement
(sp.scores.dune.ca <- data.frame(sp.scores.dune.ca))
##                   CA1         CA2
## Achimill -0.908593994  0.08460595
## Agrostol  0.933782632 -0.20651410
## Airaprae -1.004341487  3.06748567
## Alopgeni  0.400877498 -0.61838999
## Anthodor -0.966759907  1.08360766
## Bellpere -0.500177308 -0.35502842
## Bromhord -0.657624283 -0.40634288
## Chenalbu  0.424451025 -0.84401946
## Cirsarve -0.056471184 -0.76397838
## Comapalu  1.916902749  0.52150490
## Eleopalu  1.763825952  0.34562338
## Elymrepe -0.370742471 -0.74147804
## Empenigr -0.690268768  3.26420263
## Hyporadi -0.854079092  2.52821112
## Juncarti  1.275799635  0.09962851
## Juncbufo  0.081568571 -0.68074282
## Lolipere -0.502724848 -0.35955008
## Planlanc -0.840581620  0.24885595
## Poaprat  -0.389190052 -0.32998545
## Poatriv  -0.181852449 -0.53996706
## Ranuflam  1.558855994  0.30699556
## Rumeacet -0.652892711 -0.25524774
## Sagiproc  0.003639789  0.01718622
## Salirepe  0.610351084  1.54868352
## Scorautu -0.195656929  0.38883571
## Trifprat -0.881164096 -0.09792387
## Trifrepe -0.076660735 -0.02032460
## Vicilath -0.618932296  0.37139708
## Bracruta  0.182224012  0.26477423
## Callcusp  1.951985809  0.56742556

Comme on a fait acvec les scores des espèces de la PCA, on ne garde que les plus hautes valeurs absolues.

A.ca <- top_n(sp.scores.dune.ca, 3, CA1)
B.ca <- top_n(sp.scores.dune.ca, -3, CA1)
C.ca <- top_n(sp.scores.dune.ca, 3, CA2)
D.ca <- top_n(sp.scores.dune.ca, -3, CA2)
sp.scores.skim.ca <- bind_rows(A.ca,B.ca,C.ca,D.ca)
(sp.scores.skim.ca <- distinct(sp.scores.skim.ca))
##                  CA1         CA2
## Comapalu  1.91690275  0.52150490
## Eleopalu  1.76382595  0.34562338
## Callcusp  1.95198581  0.56742556
## Achimill -0.90859399  0.08460595
## Airaprae -1.00434149  3.06748567
## Anthodor -0.96675991  1.08360766
## Empenigr -0.69026877  3.26420263
## Hyporadi -0.85407909  2.52821112
## Chenalbu  0.42445102 -0.84401946
## Cirsarve -0.05647118 -0.76397838
## Elymrepe -0.37074247 -0.74147804

Graphique simple avec les fonctions plot et ordiplot.

plot(ca.dune)
ordiellipse(ca.dune, dune.env$Management, display = "sites", kind = "sd", label = T)

ggplot

On commence le graphique ggplot. On peut faire un graphique qui ressemble à celui fait par plot en ajoutant “à la main” les sites et les espèces (en fonction de leurs scores extraits).

ca.plot <- 
  ggplot(data = site.scores.ca, 
         aes(x = CA1, y = CA2)) +
  expand_limits(x=c(-4,4), y=c(-1,3)) + # sert à modifier l'échelle numérique du graphique
  geom_text(data = site.scores.dune.ca, 
            aes(x = CA1, y = CA2, label= rownames(site.scores.dune.ca)), 
            color = "black") + # ajout du texte des sites
  geom_text(data = sp.scores.dune.ca, 
            aes(x = CA1, y = CA2, label= rownames(sp.scores.dune.ca)), 
            color = "red") + # ajout du texte des espèces
  xlab("CA1 (25.34%)") + 
  ylab("CA2 (18.92%)") # contrairement avec la fonction ggord, on écrit la variation expliqué par les axes soi même
ca.plot

On peut rajouter les vecteurs de note envfit.

ca.plot <- 
  ggplot(data = site.scores.ca, 
         aes(x = CA1, y = CA2)) +
  expand_limits(x=c(-2,2.5), y=c(-1,3)) + # on zoom in sur l'échelle
  geom_text(data = site.scores.dune.ca, 
            aes(x = CA1, y = CA2, label= rownames(site.scores.dune.ca)), 
            color = "black") +
  geom_text(data = sp.scores.dune.ca, 
            aes(x = CA1, y = CA2, label= rownames(sp.scores.dune.ca)), 
            color = "red") + 
  xlab("CA1 (25.34%)") + 
  ylab("CA2 (18.92%)") +
  geom_segment(data = vect.ca.env, # ajout des vecteurs envfit
               aes(x = 0, xend = CA1, y = 0, yend = CA2),
               arrow = arrow(length = unit(0.75, "cm")), 
               color = "black") +
  geom_text_repel(data = vect.ca.env, 
                  aes(x = CA1, y = CA2, label = variable.env, size = 3), 
                  color = "firebrick1", 
                  fontface = "bold", 
                  show.legend = FALSE)
ca.plot

On va mettre des points pour les sites et changer leurs couleurs et leurs symboles en fonction de leur type d’aménagement. On va aussi mettre en graphique uniquement les espèces aux plus grands scores.

ca.plot <- 
  ggplot(data = site.scores.ca, 
         aes(x = CA1, y = CA2), 
         color = Management) + # on ajoute ici arg "color"
  expand_limits(x=c(-2,2.5), y=c(-1,3)) +
  geom_point(aes(shape = Management, color = Management), 
             size=3) + # on change de geom_text à geom_point pour les sites
  xlab("CA1 (25.34%)") + 
  ylab("CA2 (18.92%)") +
  geom_text_repel(data = sp.scores.skim.ca, 
                  aes(x = CA1, y = CA2), 
                  label=rownames(sp.scores.skim.ca), 
                  color = "black") + # on garde geom_text pour les espèces mais on change de dataframe
  geom_segment(data = vect.ca.env,
               aes(x = 0, xend = CA1, y = 0, yend = CA2),
               arrow = arrow(length = unit(0.75, "cm")), 
               color = "black") +
  geom_text_repel(data = vect.ca.env, 
                  aes(x = CA1, y = CA2, label = variable.env, size = 3), 
                  color = "firebrick1", 
                  fontface = "bold", 
                  show.legend = FALSE) +
  scale_shape_manual(values = c(15,16,17,18)) # on décide des symboles à utiliser (se référer à la charte plus haut)
ca.plot

On va ajouter des ellipses autour de nos sites en fonction de leur type d’aménagement. On va d’abord utiliser la fonction ggordiplot de la librairie ggordiplots pour extraire les données des ellipses. (on peut aussi extraire les données de hulls et spiders.) En même temps, ça nous donne un autre beau graphique, mais qui est plus complexe à modifier.

library(ggordiplots)
## Loading required package: glue
ca.ggordiplot <- gg_ordiplot(ca.dune, groups = dune.env$Management, kind = "sd", conf = 0.95, pt.size = 3) # on spécifie l'intervalle de confiance des ellipses

names(ca.ggordiplot)
## [1] "df_ord"      "df_mean.ord" "df_ellipse"  "df_hull"     "df_spiders" 
## [6] "plot"
ca.ellipses <- ca.ggordiplot$df_ellipse # extraction des valeurs des ellipses
head(ca.ellipses)
##   Group          x         y
## 1    BF -0.5850575 0.7794705
## 2    BF -0.5691512 0.7777560
## 3    BF -0.5530958 0.7726193
## 4    BF -0.5369545 0.7640805
## 5    BF -0.5207911 0.7521735
## 6    BF -0.5046694 0.7369452

On prend les valeurs des ellipses qu’on vient d’extraire pour mettre en graphique dans ggplot ces ellipses. Om change aussi les couleurs de nos groupes.

ca.plot <- 
  ggplot(data = site.scores.ca, 
         aes(x = CA1, y = CA2), 
         color = Management) + 
  expand_limits(x=c(-2,2.5), y=c(-1,3)) +
  geom_point(aes(shape = Management, color = Management), 
             size=3) +
  xlab("CA1 (25.34%)") + 
  ylab("CA2 (18.92%)") +
  geom_text_repel(data = sp.scores.skim.ca, 
                  aes(x = CA1, y = CA2), 
                  label=rownames(sp.scores.skim.ca), 
                  color = "black") +
  geom_segment(data = vect.ca.env,
               aes(x = 0, xend = CA1, y = 0, yend = CA2),
               arrow = arrow(length = unit(0.75, "cm")),
               color = "black") +
  geom_text_repel(data = vect.ca.env, 
                  aes(x = CA1, y = CA2, label = variable.env, size = 3), 
                  color = "firebrick1", 
                  fontface = "bold", 
                  show.legend = FALSE) +
  scale_shape_manual(values = c(15,16,17,18)) +
  geom_path(data = ca.ellipses, 
            aes(x = x, y = y, color = Group), 
            show.legend = FALSE) + # ajout des ellipses
  scale_color_brewer(palette="Dark2") # changement des couleurs des ellipses et des sites
ca.plot

Enfin, on pourrait rajouter notre thème personnalisé, mais optons pour un existant comme “bw”.

ca.plot + theme_bw()

Comme mentionnée précédemment, vous pouvez modifier presque à l’infini votre graphique!

On sauvegarde notre figure de CA.

ca.plot + theme_bw()

ggsave("ca.png", path="./figures/")
## Saving 7 x 5 in image

PCoA

Préparation des données

Faire une PCoA avec la fonction prcomp.

dune.bray <- vegdist(dune, method='bray')

Faire une PCoA dans R

pcoa.dune <- cmdscale(dune.bray, k =(nrow(dune) - 1), eig = TRUE)
## Warning in cmdscale(dune.bray, k = (nrow(dune) - 1), eig = TRUE): only 14 of the
## first 19 eigenvalues are > 0
summary(pcoa.dune)
##        Length Class  Mode   
## points 280    -none- numeric
## eig     20    -none- numeric
## x        0    -none- NULL   
## ac       1    -none- numeric
## GOF      2    -none- numeric

Dans cet exemple, on va s’intéresser aux axes 1 et 3, pour faire changement!

Extraire les proportions expliquées par les axes d’intérêt (1 et 3).

pcoa.eig2 <- eigenvals(pcoa.dune) # extraction des eigenvalues
pcoa.eig2 <- abs(pcoa.eig2) # calcul de leurs valeurs absolues
sum.pcoa.eig2 <- sum(pcoa.eig2) # total des eigenvalues
pcoa.prop.exp <- pcoa.eig2/sum.pcoa.eig2 # proportion expliquée par chacun des axes
(pcoa.axis1 <- pcoa.prop.exp[1])
## [1] 0.3510557
(pcoa.axis3 <- pcoa.prop.exp[3])
## [1] 0.09439071

Extraire les scores de la PCoA.

site.scores.pcoa <- scores(pcoa.dune) # site scores
sp.scores.pcoa <- wascores(pcoa.dune$points, dune) # species scores 
site.scores.pcoa <- data.frame(site.scores.pcoa)
sp.scores.pcoa <- data.frame(sp.scores.pcoa)
Management <- dune.env$Management
site.scores.pcoa <- cbind(site.scores.pcoa, Management)

Comme on a fait acec les scores des espèces de la PCA et de la CA, on ne garde que les plus hautes valeurs absolues, mais cette fois ci, des axes 1 et 3!.

A.pcoa <- top_n(sp.scores.pcoa, 3, X1)
B.pcoa <- top_n(sp.scores.pcoa, -3, X1)
C.pcoa<- top_n(sp.scores.pcoa, 3, X3)
D.pcoa <- top_n(sp.scores.pcoa, -3, X3)
sp.scores.skim.pcoa <- bind_rows(A.pcoa,B.pcoa,C.pcoa,D.pcoa)
(sp.scores.skim.pcoa <- distinct(sp.scores.skim.pcoa))
##                    X1            X2           X3          X4           X5
## Comapalu  0.458911230  0.1005585008  0.125102902 -0.10064941 -0.116930204
## Eleopalu  0.439782223 -0.0007156372  0.106965962 -0.02841692 -0.008577587
## Callcusp  0.478165876  0.0484426585  0.130601445 -0.08277595 -0.033978450
## Achimill -0.288291335  0.0379745204  0.024675952 -0.09206847 -0.054125785
## Bromhord -0.257640216 -0.0863954310  0.007271944 -0.04584125 -0.050452440
## Trifprat -0.277616388  0.0548361248 -0.015862669 -0.04290624 -0.172247260
## Chenalbu  0.179872241 -0.2151706072 -0.273949689 -0.11812429  0.054398976
## Empenigr -0.005394285  0.4691160578 -0.217330420  0.14961117  0.161953074
## Juncbufo  0.059380874 -0.1672334185 -0.216754109  0.01068520 -0.006485779
##                   X6           X7           X8           X9          X10
## Comapalu -0.08308461 -0.034691754 -0.006325968  0.046184733 -0.002246251
## Eleopalu  0.03187448  0.003764263  0.018051782 -0.005171053 -0.011032580
## Callcusp -0.02214687  0.023556133  0.013738043 -0.002484335  0.032798306
## Achimill -0.00227727  0.014193770  0.015356868  0.001709240 -0.008461562
## Bromhord -0.05616890  0.070270108  0.040802596  0.013692138 -0.034231159
## Trifprat  0.15050556 -0.017699829  0.014358724 -0.002775301  0.025056357
## Chenalbu -0.09378649 -0.035811438 -0.138395894 -0.091709220 -0.037094509
## Empenigr -0.09775372  0.044962813  0.071770276  0.049827956  0.025898170
## Juncbufo  0.03276354 -0.054325318 -0.034914441  0.026597573  0.011229702
##                   X11          X12           X13          X14
## Comapalu  0.072444555 -0.006104436  0.0127702458  0.007221267
## Eleopalu -0.008371554 -0.014122056  0.0008705509 -0.002602251
## Callcusp -0.035297813  0.016298900 -0.0011275608  0.002066468
## Achimill -0.012001757 -0.001817249 -0.0039917245 -0.004955411
## Bromhord -0.007920231  0.005418645  0.0040940773 -0.003890661
## Trifprat -0.028308593 -0.011715116  0.0154123351  0.010475144
## Chenalbu -0.040870761 -0.002477992  0.0321341180  0.003153000
## Empenigr -0.041757420 -0.042319419 -0.0135562316  0.017897602
## Juncbufo -0.004033853  0.018627584  0.0021911465 -0.006478585

envfit SUR LES AXES 1 ET 3 (on n’oublie pas!)

(pcoa.envfit <- envfit(pcoa.dune, dune.env, choices=c(1,3)))
## 
## ***VECTORS
## 
##         Dim1      Dim3    r2 Pr(>r)  
## A1 0.9999700 0.0078085 0.379  0.018 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
## 
## ***FACTORS:
## 
## Centroids:
##                 Dim1    Dim3
## Moisture1    -0.2650  0.0607
## Moisture2    -0.1620  0.0000
## Moisture4     0.1065 -0.2369
## Moisture5     0.3271  0.0069
## ManagementBF -0.2694  0.0706
## ManagementHF -0.1350 -0.0394
## ManagementNM  0.1822  0.0370
## ManagementSF  0.0650 -0.0395
## UseHayfield  -0.0608 -0.0240
## UseHaypastu  -0.0227 -0.0210
## UsePasture    0.1213  0.0673
## Manure0       0.1822  0.0370
## Manure1      -0.1654  0.0175
## Manure2      -0.1648 -0.0944
## Manure3       0.1397 -0.0451
## Manure4      -0.1656  0.0944
## 
## Goodness of fit:
##                r2 Pr(>r)    
## Moisture   0.6918  0.001 ***
## Management 0.2634  0.130    
## Use        0.0614  0.669    
## Manure     0.2892  0.196    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999

Pour incorporer les vecteurs résultant de ce test, il faut d’abord les extraire et les mettre dans un dataframe.

vect.pcoa.env <- data.frame(pcoa.envfit[["vectors"]][["arrows"]])
vect.pcoa.env$variable.env <- rownames(vect.pcoa.env)
vect
##          PC1       PC2 variable.env
## A1 0.9722159 0.2340862           A1

PERMANOVA. On va tester si la composition végétale diffère entre différents type d’aménagement (facteur; variable quantitative) des sites.

dune.dist.bray <- vegdist(dune, method='bray')
disper.dune.bray <- betadisper(dune.dist.bray, dune.env$Management)
anova(disper.dune.bray)
## Analysis of Variance Table
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq F value Pr(>F)
## Groups     3 0.13831 0.046104  1.9506 0.1622
## Residuals 16 0.37816 0.023635
permanova.dune.pcoa <- adonis2(dune.dist.bray ~ Management, data = dune.env)
permanova.dune.pcoa
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = dune.dist.bray ~ Management, data = dune.env)
##            Df SumOfSqs      R2      F Pr(>F)   
## Management  3   1.4686 0.34161 2.7672  0.004 **
## Residual   16   2.8304 0.65839                 
## Total      19   4.2990 1.00000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Extraction des ellipses de la PCoA (axes 1 et 3!)

library(ggordiplots)
pcoa.ggordiplot <- gg_ordiplot(pcoa.dune, choices = c(1,3), groups = dune.env$Management, kind = "sd", conf = 0.95, pt.size = 3)

names(pcoa.ggordiplot)
## [1] "df_ord"      "df_mean.ord" "df_ellipse"  "df_hull"     "df_spiders" 
## [6] "plot"
pcoa.ellipses <- pcoa.ggordiplot$df_ellipse

ggplot

pcoa.plot <- 
  ggplot(data = site.scores.pcoa, 
         aes(x = Dim1, y = Dim3), 
         color = Management) + 
  expand_limits(x=c(-0.35,0.6), y=c(-0.35,0.35)) +
  geom_point(aes(shape = Management, color = Management), 
             size=3) +
  xlab("PCoA1 (35.10%)") + 
  ylab("PCoA3 (9.44%)") +
  geom_text_repel(data = sp.scores.skim.pcoa, 
                  aes(x = X1, y = X3), 
                  label=rownames(sp.scores.skim.pcoa), 
                  color = "black") +
  geom_segment(data = vect.pcoa.env, 
               aes(x = 0, xend = Dim1, y = 0, yend = Dim3),
               arrow = arrow(length = unit(0.75, "cm")), 
               color = "black") +
  geom_text_repel(data = vect.pcoa.env, 
                  aes(x = Dim1, y = Dim3, label = variable.env, size = 3), 
                  color = "firebrick1", 
                  fontface = "bold", 
                  show.legend = FALSE) +
  scale_shape_manual(values = c(15,16,17,18)) +
  geom_path(data = pcoa.ellipses, 
            aes(x = x, y = y, color = Group), 
            show.legend = FALSE) + 
  scale_color_brewer(palette="Dark2") +
  super_theme
pcoa.plot

ggsave("pcoa.png", path="./figures/")
## Saving 7 x 5 in image